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1. Introduction. In his recent work in Generalized Axially Symmetric Potential 
Theory, A. Weinstein [1] has pointed out that the flow about an axially symmetric body 
in ordinary three-dimensional space may be obtained from the electrostatic potential of 
the five-dimensional body of revolution possessing the same meridian profile. This 
method of solution is referred to as the method of Generalized Electrostatics. It has been 


recently used by L. E. Payne and A. Weinstein [2] in deriving a relationship between 


capacity and virtual mass and has been employed by A. Weinstein [3] in solving certain 


torsion problems. 

In this paper the method of Generalized Electrostatics is used to obtain the flow 
about a spindle and a lens. The flow about a spindle seems not to have been treated in 
the literature despite its obvious importance. The lens problem however was treated in 
1947 by M. Shiffman and D. C. Spencer [4] who applied an ingenious and difficult pro- 
cedure involving the method of images in a multi-sheeted Riemann-Sommerfeld space. 
The solution given in this paper is a straightforward generalization of results given already 
in 1868 by F. G. Mehler [5]. Our formulas are considerably simpler than those of Shiffman 
and Spencer but there is no obvious computational method of showing that their solution 
may be reduced to ours or vice versa. However, the identity of these two solutions is 
guaranteed by a uniqueness theorem. In an oral communique Professor D. C. Spencer has 
pointed out the fact that the problem of the spindle would be difficult if not impossible to 
solve by the method of images in a Riemann-Sommerfeld space. It will be seen that 
generally speaking little difficulty is encountered in extending to an arbitrary odd-dimen- 
sional space the known solutions of three-dimensional electrostatics problems. Hence by 
the method of Generalized Electrostatics we obtain the flows about axially symmetric 
profiles almost immediately from the electrostatic solutions. 

In this paper we are concerned chiefly with bodies of revolution in three- and five- 
dimensional spaces; but since spaces of other dimensions have various applications we 
shall first obtain the solution in a general n-dimensional space. The ordinary three-dimen- 
sional terminology will be retained throughout. Later we shall assign the particular value 
to n which is demanded by the physical problem. 

2. Basic Equations. We shall restrict ourselves in this paper to profiles of revolution 
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in a uniform stream. We shall assume further that the fluid is incompressible and that at 
infinity it is travelling at uniform velocity U parallel to the axis of symmetry. Let xy be 
the meridian plane of an n-dimensional body which is symmetric about the z-axis. We 


shall define all functions in the meridian half plane y = 0. 
An axially symmetric potential function g{n} in a space of n-dimensions is a solution 


of the partial differential equation 


a 2 oy{n} ra) doin} E 
- ly et) + - (y a, Q, n & 2. (1) 
, Ox oy \’ OY 


The associated stream function ~{n} is defined with the aid of the generalized Stokes- 


Beltrami equations 


Opin} Opin} n-2 O9{n} dyin} (eo 
i > . ’ y : _ = ’ (4) 
Ox OY : OY Ox 
Let V{n} Uy" '(n 1)"" — W{n} be the stream function describing our flow. We 
assume V {7} to vanish on the profile and along the z-axis.* We make use of a correspon- 
dance relationship [1], namely y{n} Uy""(n — 1)"'e{n + 2}, to obtain the funda- 
mental equation 
Vin} Uy" '(n 1)-'(1 — of{n + 2}). (3) 


This equation relates the stream function V{n} for an n-dimensional body of revolution B 
+ 2} in a space of n + 2 dimensions. The 


to an electrostatic potential function ¢{n 
potential g{n + 2} assumes the value unity on the profile boundary and vanishes at 
infinity. 


We may by introducing the substitution x{n} y'" ”““of{n} reduce Eq. (1) to the 


form 

Vx [(n 2)(n — 4)/4y" |x 0 (4) 
where VY’ denotes the Laplacian operator. It is obvious that for n = 2 or n 1, x is 
harmonic. C. Snow [6] in a discussion of non-axially symmetric potential problems in 
three dimensions has treated solutions of Eq. (4) for odd values of n. 
f(& + in) Eq. (4) takes the form 


Under a transformation x + iy = z = f(£) 
Kee + Xe [(n 2)(n — 4)/4h*y*]}x = 0 (5) 
where h’ f’(¢) |~*. Similarly under such a transformation Eq. (1) becomes 
y ; 1 4 f i : 
- (y: ° 4 | = (y" = ae) = (). (6) 
0& \' 0& On \’ on 
If y is of the form f(t) - g(m) the solution of (6) is readily obtained by separation of 
variables. If on the other hand (h’y’)~’ f,(€) + g,(m) as in the cases considered here, 


we find it more convenient to use Eq. (5). 
Payne and Weinstein [2] have derived a relationship between the n-dimensional 
virtual mass M{n} (uniform flow in the x-direction) and the n + 2-dimensional capacity 
{n + 2} which for n 3 is given by 
M{3} + V{3} = (2r/3)C{5}. (7) 


*This implies that the profile and the z-axis form a single streamline. 
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In (7) V{3} denotes the volume of the body of revolution (the density of the fluid has 
been taken as unity). The capacity C{n} is determined from the equation 
Cin} (n 2) lim r"“g{n}, rP=o2t+y’ (8) 
where g{n} is the electrostatie potential of the body. We shall make use of (7) and (8) 
in determining the virtual mass of the spindle and the lens. 
3. Flow About a Spindle. A curve = &, in dipolar coordinates defines the profile of a 
a spindle. The dipolar transformation is given by 


x+y ic cot 3(E + tn) (9) 
where ¢ is a positive constant. The range of coordinates is chosen as — © < » <+®, 
0) é < wr. The boundary of the spindle is given by — = & < m and the exterior region 
is defined by 0 < & < & . We may choose as particular solutions of Iq. (1) functions 
determined with the aid of (5) which are of the form 

l(s — B/A (yi ye'2ai™ (id) cos my (10) 

where 2¢ ! 3, cosh n, 4 = cos — and X = cot & The Q function represents a 

generalized spherical harmonic of the second kind as defined by I. Hobson [7]. The fune- 

tions considered in (10) obviously vanish at infinity (& = 0, » = 0). We are primarily 

interested in odd integral values of n and in this case (10) is considerably simplified. 
Particulir solutions of (1) are then given by 

(s t)""'? KS? (t) cos mn (11) 


where (g) denotes the gth partial derivative with respect to the argument and K,(0) is a 


Legendre funetion of complex degree commonly called a conal function [7, p. 445]. It is 


defin al 
2 J 1/2 
KC ( } cosh amr | {2 cosh wu + 2 cos é cos au du. (12) 
\T 0 
If we replace — by (# — &) in (12) we obtain 
/9 ; a 
K i) = cosh am | (2 cosh u — 2 cosé]’~ cos au du. (13) 
\r Jo 
Now if 0 < & < x, Eq. (13) may be differentiated q times with respect to cos — the order 


of differentiation and integration being interchanged on the right. For 0 < & < m the 
term (cosh » — cos £9) satisfies the conditions of the Fourier integral theorem and 


may be expanded as follows 


(cosh ” COS & 


a : 
(3) | cos anh | [cosh n’ — cos &] ‘7"'’” cos an’ in’) da. (14) 


We notice that the term in braces is simply K"’(—2) multiplied by a function of a. We 
have used f) to represent cos £ . We understand by the superscript (q) in K{(—d) the 
gth partial derivative with respect to cos — and not with respect to —cos &. 

We now choose for the electrostatic potential g{n} a function of the form 


a 


gin\(s — 2°" = | A(a)K ,"'(t) cos an da. (15) 


“0 
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| for € = & determines with the aid of (14) the function A (a) 


So 


The condition that ¢ 
since (15) must be satisfied identically in ». The n-dimensional electrostatic potential 
g{n} is found for odd n to be 

© Kh (— bh ) Ke (2 eos ax . 
| - I da. (16) 
“0 K'(t.) cosh ar 


€ 


This integral converges for all 7 and for 0 S$ — < & 
1 odd integer the electrostatic potential function becomes much more 


If n is not ar 


complicated in general. However, in the special case n = 4 the potential may be easily 
determined with the aid of (5) as 
, f° sinh a(r — &) sinh at cos an = 
of4} = %s— Na — ef)? | a de. (17) 
sinh az sinh a&, 


For odd values of n the stream function V{n} representing the flow about a spindle may 


now be obtained from (16) with the aid of (3). We have 


U(e sin & 
Vin} = = 
Aq + 1\(s t 
[ (Qe)'/*(g — #)° [ K'""'(—t) K5""" (8) cos an ] 18 
: aa : - aa |. \ ) 
0 l(q + 3/2) ‘ ty "(¢t,) cosh am al 
Thus when n = 3 (g = 0) we obtain the stream function corresponding to a uniform flow 
about a three-dimensional spindle. It is given by 
Uc" sin’ é : os Me cet Kh eo a 
¥i3} = > >} 1 — (2) "(s — by” | ear - - + ° (19) 
a\5 t “0 K.’(t)) cosh ar 
If & = 2/2 the spindle becomes a sphere and the well known stream function for the 


sphere is obtained. 
The capacity C{n} of an n-dimensional spindle is found according to (8) as 
) 9 
da. (20) 


[ Ka (—t)Ka (1 


K ."'(t,) cosh ar 


The case n = 3 has been given by G. Szegé [8]. We obtain the virtual mass 4/ {3} with 


the aid of (7 


P 
/ <>’ 
M{3} = —(<)mc\2 + 3 cot? & + 3 cot & esc? é 
/ 
(da? + IKI (—t) , | 
+3 | ~ da?. (21) 
Jo K.. (to) cosh ar ) 
It is easily verified that for & = 2/2 we have the well known virtual mass of the sphere. 


4. Flow About a Lens. Let us introduce the peripolar transformation, x + 7y 
1(¢ + in), where c is a positive constant. The profile of a lens is defined by two 


—C Col 
We shall assume that 0 < &, < & < 2x. The external region 


curves & = £, , and é& = &,. 
is chosen as n > 0, & — 2m < & < &,. Particular solutions of (1) are given by functions 


of the type 


(gs — t)**""*(s’ — 1)°”°K 3(s) cosh mé (22) 
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and those obtained by replacing cosh mé by sinh mé. In (22) KA(s) is a Legendre function 
of the type discussed by Mehler, see [7, p. 451]. As in the case of the spindle the lens 
problem is solved much more easily and the solution is given in a much simpler form 
whenever v is an odd integer. Particular solutions of (1) are in this case given by 


(s — t)**' K‘2(s) cosh mé. (23) 
solutions containing sinh mé being understood as before. The function K,(s) is defined 
fs p. 151] as 


" - 
K..(s) = (2) cosh amr | [2 cosh u + 2 cosh n]'”’ cos au du. (24) 


T “0 


We have aiso the known expansion [7, pp. 452, 453] valid for 0 < —& < 2x 
@-)'? =?” cosh a(é — m) sech am K,(s) da. (25) 


Clearly K .’’(s) is a well defined function obtained from (24) by a permissible exchange 
of order of integration and differentiation. Also for 0 < — < 2x Eq. (25) may be differen- 
tiated g times with respect to s the order of integration and differentiation being enter- 
changed on the right. We choose as the electrostatic potential g{n} of the lens 


| 


gini(s — t) “ a a | [A(a) cosh af + B(e) sinh aé|K;," (s) da. (26) 
The functions A(a@) and B(a) may be chosen in such a way that g{n} = 1 foré = é, and 
£ = §. We need only differentiate (25) ¢ times with respect to s evaluate for § = £, and 


£ = £, and insert in (26). Since (26) must be satisfied identically in y the functions A (a) 
and B(a) are easily determined. The electrostatic potential of an n-dimensional lens 
(n odd) is thus given by 
gin} = (—1)%(2n)'?T '(q + 1/2)(s — 7°"? | F(a, ¢) sech aw K‘?’(s) da (27) 
“0 


W here 


sinh a(2r — &, 4. £,) F(a, &) 
sinh a(é, — £) cosh a(x — &) + cosh a(r — ~,) sinh a(2r — & + &). (28) 


Equation (27) is valid for all positive 7 and for all in the interval £, — 27 < — < &, . The 


case n = 3 was given by F. G. Mehler [5] in 1868. 


We note again that in case n = 4 the electrostatic potential may be easily obtained 
with the aid of (5). We have 


~~ 


{4} = 2(s — d/(s’ — 1)'”] | 


“0 


F(a, &) esch az sin an da. (29) 
The stream function ¥{n} representing the flow about an odd-dimensional lens is 


obtained from (27) with the aid of (3). Thus for n = 3 


W{3} = [Uc'(s’ — 1)/2(s — t)’] E + 2°75 — t)*” [ F(a, ¢) sech aw K{’(s) ia, (30) 
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By an appropriate choice of £, and & electrostatic potentials and stream functions for a 
spherical bowl, symmetrical lens, hemisphere, e/c., may be determined. 


Jecause of the invariance of form of the Stokes-Beltrami equations, 1.e. 
y © = V,, yn @. WV 31) 


it isa simple matter to determine from ¥ the velocity potential & to which ¥ is associated: 
It remains to be shown that this soluion is unique. 

The problem of establishing the uniqueness of the stream function Y defined in the 
infinite region outside the profile in the zy plane and satisfying prescribed boundary 
conditions is equivalent to the problem of establishing uniqueness of this function in an 
infinite half strip in the én plane. An application of Green’s formula would demand a 
knowledge of the behavior of the derivatives of VY at infinity in the & plane. Hence we 
find it more convenient to make use of an eigen value method due to A. Weinstein [9] 
which requires only a knowledge of the behavior of V at infinity. By this method we can 
show that there is only one stream function V which satisfies the prescribed conditions on 
the lens profile. If the stream function is unique the potential ® is also unique up to an 
arbitrary constant which must be zero in order that the potential vanish at infinity. 

The electrostatic capacity of an n-dimensional lens is determined for integral values 
of q (n odd) with the aid of Iq. (8). It is given by 


Z l/&@ 2g I }¢ 
ear | 
Cyn : 
l'(q L/2 
[ sinh ake, ( ish a\T So) + cosh CAT — &,) sinh a(2r - £,) ite ] 99 
| ‘ i YS (1) da. OZ) 
sinh a(2xr + ¢ £) cosh am 
The case n = 3 has been given by G. Szego [8]. The virtual mass 1/{3} is obtained with 
the aid of [7] and given by 
ria) 9 *® sinh ag, cosh a(r &) + cosh alr — &,) sinh a(2r — é 
M {3} Qa : ; 
sinh a(2r + é, £,) cosh amr 
(4a° + 1) da — V}3!} 33) 
where 
V $3} re /6)4(2 cos €,) cot —,/2 ese” &,/2 — (2 cos £) cot &/2 ese” &/2 


Equation (33) is a much simpler expression for the virtual mass than that given by 
Shiffman and Spencer. Several special cases may be obtained easily from (33). In particu- 


lar the virtual mass of a hemisphere is given by 
M3! = (2mc' /81)(135 59(3) "a = 2.545¢ 34) 


where c is the radius of the hemisphere. 

5. Additional Results. We shall list here the electrostatic potentials of certain other 
n-dimensional bodies of revolution. With the aid of (3), (7) and (8) the corresponding 
flow problems can be completely solved. The results given in this section are valid for 
any positive real value of g(n > 2). It will be noted that the results are simplified con- 


siderably whenever n is an odd integer. 
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I. Sphere: The potential of an n-dimensional sphere of radius a is given as 


g{n} = (a/r)"~ r=(2° + y)'”. (35) 


’ 


Il. Two separated spheres: The lines » = a < O and ny = 6 > 0 in dipolar coordinates 
define two separated spheres. The potential in this case is given by 


NB . r Na . , 
: 7 Tf sinh N(n —a)t+e sinh N(8 — n) P.(t/2q + 1) (36 
rr ™ a sinh N(gB — a) i d } 


where V = n+94-+ and P(t | 2q¢ + 1) represents the 2g + 1-dimensional zonal spherical 
harmonic or more commonly the Gegenbauer polynomial. 
Ill. Prolate Spheroid: A line — = & defines a prolate spheroid under the transformation 
c cosh ¢. The electrostatic potential for such a spheroid is 


ein} = (po/p)"Q%s)/Q2(s) (37) 
where p = sinh mm 4 cosh oe 
IV. Oblate Spheroid: Under the transformation z = c sinh ¢ an oblate spheroid is 
defined by a line é £,, and the potential is given as: 
gin} = (S85 8)“Qi(ip) ()"(1po). (38) 
V. Dise: lf i 0 the oblate spheroid becomes a disc of radius c and the potential of 


such a disc is obtained from (38) as 


yin} | 2 ' exp {(2+ ‘ei r(« + Yo(d)e QQ (tp). (39) 


We have listed here only a few examples. Numerous others can be easily obtained. 
6. Internal Pioblems. The method of Generalized Electrostatics is also useful in 
determining the flow induced in a fluid contained between two or more axially symmetric 
boundaries when one or more of the boundaries moves with respect to the others at uni- 
form velocity parallel to the axis of symmetry. In this case we consider instead of Eq. (3) 


the equation 


i , , 

Vin} = Uy" '(n — 1) 'ef{n + 2}. (40) 

On a moving boundary ¥{n} = V,y""'(m — 1)7' where V; = ¢,U (ce; is a constant possibly 
differing for each moving boundary). On a stationary boundary ¥{n} = 0. This problem 


is reduced by (40) to the solution of a steady state heat flow problem in n + 2-dimensions. 
The boundaries in n + 2-space corresponding to the moving boundaries in n-space are 
maintained at temperatures c, and those corresponding to the stationary boundaries are 
kept at temperature 0. This procedure applies in particular to the case in which the fluid 
is bounded by two eccentric spheres, two tori or to the ease in which one portion of the 
boundary moves with respect to another portion in an infinite fluid. 

7. Concluding Remarks. In this paper we have considered only three dimensional 
flow problems. It should be remarked, however, that a similar procedure may be em- 
ployed in solving plane flow problems for profiles symmetric about the z-axis. In fact if 
the profile possesses symmetry with respect to both axes the plane problem may be 


solved for uniform flow in any direction. 
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A FORMULA FOR AN INTEGRAL OCCURRING IN THE THEORY OF 
LINEAR SERVOMECHANISMS AND CONTROL-SYSTEMS* 


BY 
HANS BUCKNER 


Minden, Germany 


Introduction. Let ¢ denote the time, p = d/dt the differential operator with respect 


to time and 
f, (vt) = ax" + av” +--+ +4, ; a ~0,n>1 (1) 
a polynomial with real coefficients. If all zeros of f,(x) have negative real parts, every 


solution y(t) of 
f.(p)y = 0 (2) 


and all derivatives py tend to zero with increasing ¢. Moreover the integral 


Y= | y(t) dt (3) 
exists. The purpose of this paper is to develop a formula for Y in terms of squared linear 


forms of the initial values 
p yO) = ” 3 k=Q@Q,1, ++ .n— 1. (4) 


No further quantities but the coefficients a; of (1) shall appear in this formula. 
Such a formula may be useful for the design of linear servomechanisms and control- 


systems, governed by the equation 
f.A(p)y = 2(0). (2’) 


where z(f) may be considered as an arbitrary disturbance function. For instance, let 
z(t) 1 for < 0. At ¢ = 0, z(t) may step down to z(t) = 0 for ¢ > 0. The response 
y(t) then is a solution of (2), and the integral Y measures, how fast the systems lines 
up with the stepping of z. The knowledge of Y makes it possible to choose the coefficients 
a, of (1) under given conditions in order to minimize Y**. Two examples of such a mini- 
mization will be given in Sec. 4. 

The development of this formula will also yield a new approach to the well known 
Hurwitz criterion of stability and to reductions of ‘stable’ operator polynomials in p 
to such of a lower degree, including the reduction of Schur [1]. 

1. Auxiliary theorems and algorithms of reduction. Notation. Let J be the imaginary 
axis of the complex plane, J’ the set of all points wi, of J with w > 0, J” the set of all 
points wi of J with w < 0, and Re z the real, Im z the imaginary part of z. 

Definitions. Let f(x) = box” + ba" ' + +++ + b,, a polynomial with real or complex 
coefficients. We call m the proper degree of f(x), if b) # 0. We define now 
1) F(x) = box” + b,, as the “simplification” of f(x), if f(z) has the proper degree m, 


*Received August 3, 1951. 
**Minimization of Y has already been investigated by P. Hazebroek and B. L. van der Waerden [2] 
who also gave a formula expressing Y as a symmetric function of the zeros of (1) for special systems (4). 
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2) g(x) = b, + Dn-2x” + --- as the even and h(x) = f(x) — g(x) as the odd component 
of f(x 

3) f(x) as a “Hurwitz-polynomial’’, if all zeros of f(x) are in the left-hand half-plane 
Re x < 0 (the case f(x) = const. ~0 to be included), 


4) f(x) as definite (semidefinite) on a given set M of points of the complex plane, if a 
suitable constant c ~ 0 can be found, so that cf(x) > 0 (20) on M (for instance, 


x” is definite on J’), c to be normalized by | c| = 1. 
Lemma 1. Let p(x) and q(x) be two polynomials. The linear combination r(x, f) = 
tp(x) + (1 — #)q(x) shall have proper degree m for all values 0 < ¢ < 1. We further 


assume 7(z, f) * 0 on J for all these values of ¢. Then p(2) and q(x) have the same 
number of zeros for Re x > 0 and for Rez < 0. 

Proof. No zero of r(z, t) can pass J or can go to infinity, when ¢ is running from 0 
to 1. Hence the number of zeros for Re x > 0 remains constant. The same holds for 
Rez < 0 

Lemma 2. Let f(x) = box” + +--+ + b,, be a polynomial with real coefficients b, 
The proper degree shall be m > 1; f(x) and its simplification F(x) shall not vanish on 
J. Then f(x) and F(z) have the same number of zeros for Re x > O and for Re x < 0, 
if at least one of the following conditions is satisfied: 

a) the even component g(x) of f(x) is semidefinite on J’; 
b) m is odd, and the odd component h(x) of f(x) is semidefinite on J’; 
c) mis even, and the odd component h(x) of f(x) is definite on J’. 

Proof. We set r(z, t) = tf(z4) + 1-2 F(a) = bor” + tb,x” ; +--- +b, ,7~+b,,. 
The proper degree of r(x, ¢) is m for ali values of ¢. We shall prove that r(x, ) # 0 on 
J for 0 < ¢ < 1. Application of the first lemma then completes the proof. 

From the assumptions it follows that r(z, 0) # 0 on J and r(z, 1) ¥ 0 on J; further- 


more 7r(0, t) = b,, + 0. It is therefore sufficient to prove r(x, t) # 0 on J’ or J” for 
0 < t < 1. We denote by G(x) the even and by H(z) the odd component of F(x). G(x) 
is either the simplification of g(x) or equal to g(0) = b,, ; H(z) is either the simplification 
of h(x) or equal to h(0) = 0. 

If any polynomial s(z) is semidefinite on J’, we have cs(x) > 0 on J’ with a suitable 
constant c (!¢| = 1). Considering extremely small and extremely great values of | x 


we find cS(x) > 0 on J’ with the same constant c for the simplification S(x) of s(x). 


With this in mind, we distinguish the following three cases according to the conditions 


a, b, c of the lemma. 


a) g(x) is semidefinite on J’. This leads to cg(x) > 0 and to cG(x) > 0 on J’. We have 
either G(x) = b,, or G(x) = F(x), and in both cases G(x) ¥ 0 on J. Therefore, 
r(x, t > | Rer(z, t) | = c.tg(xz) + ¢.(1 — )G(x) > (1 — YcG(x) > Oon J’. 

b) m odd, h(x) semidefinite on J’. We have H(xr) = boa” ¥ 0 on J’ and a suitable 


constant c, making ch(x) > 0 and cH(x) > 0 on J’. Hence for 0 < t < 1 on J’ 


r(x, t)| > | Im r(a, t) | = eth(z) + ¢.1 — )H(x) > c(1l — OH(z) > 0. 


“ 


c) m even, h(x) definite on J’. We find | r(z, t) | > t| h(x) | > Oon J’ forO < ti <1. 
Thus, r(x, t) # Oon J. 

Lemma 3. Let p(x) and q(x) be any two polynomials with real coefficients, p having 
proper degree m and g having proper degree m’ < m. The polynomial f(x) = p(x)q(—2) 
and its simplification F(x) shall not vanish on J; f(x) and F(x) shall have the same 
number of zeros for Re x > 0 and for Re x < 0. From this it follows that 
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a) if p(x) is a Hurwitz-polynomial, g(x) is also one with m’ = m — 1, 
b) if g(x) is a Hurwitz-polynomial, if furthermore m = m’ + 1, and if all coefficients 


of p(x) are positive, then p(x) is also a Hurwitz-polynomial. 
Proof. The number of zeros of F(x) in the half-plane Re x < 0 may be n, the number 


of zeros in Re x > 0 may be n’. All zeros of F(x) form a regular polygon for n + n’ > 3, 


and no zero can appear on J. Hence |n — n’| < 1. Should p(x) be a Hurwitz poly- 
nomial, f(x) and F(x) have at least m zeros in Rex < 0 and not more than m’ < m — 1 
zeros in Re x > 0. Therefore n = mand n’ = m’ = m — 1. The m — 1 zeros of f(z) 
in Re x > O are those of g(—x). This means, that g(x) is a Hurwitz polynomial. Should 
the conditions of b) be satisfied, then at least m — 1 zeros of F(x) and consequently 
of f(x) appear in Re x < 0. Therefore p(x) has m — 1 zeros in Re x < 0. Should the 


last zero of p(x) be situated in Re x > 0, it must necessarily be real, i.e. positive. But 
no such zero can exist, since p(x) is assumed to have positive coefficients. This completes 
the proof of the lemma. 

Note. The condition, that all coefficients of p(x) are positive is—apart from a constant 
factor—a necessary condition for p(x) to be a Hurwitz polynomial. It is well known and 
it ean easily be proved by splitting p(x) into root factors. No coefficient can vanish 
without reducing the degree of the polynomial. 

Algorithms can be based on Lemmas 2 and 3 in order to reduce a Hurwitz polynomial 
to such of a lower degree. It may be worthwhile to explain, how the well-known reduction 
of Schur (see [1]) can be obtained in this way. 

Schur’s algorithm of reduction. We consider the polynomial (1) with real coefficients, 
but we do not assume that it is a Hurwitz polynomial. We denote by g* the even and 
by h* the odd component of f,,(x). With Schur we introduce 


fn—-(x) = (2a, — aox)[g* (x) + h*(x)] + (—1)"aor[g"(x) — h(a] (5) 
with lower degree than n. The odd component of the polynomial f(x) = f,(x)f,-1(—2) is 
h(x) = —2a,xh**(x) for even n, h(x) = 2a xg™*(x) for odd n. (6) 


This component is obviously semidefinite on J’ and on J”. It can easily be seen, that 


f,(2) = 0on J in any point z leads to f,_,(x) = 0 for the same point. Vice versa, a,f,(x) = 
0 is a consequence of f,-,(2) = 0 in any point x of J. We now assume that 
a, # 0. (7) 
This is necessary and sufficient for f,-, to have the proper degree n — 1. The polynomial 
f(x) then has the proper degree 2n — 1. If either f, or f,-; is a Hurwitz polynomial, f 
cannot vanish on J. Also F(x), the simplification of f, cannot vanish on J. Hence Lemma 
2 is applicable to f and F and the Lemma 3 to f, and f,-, . Thus, if f, is a Hurwitz 
polynomial, f,-, is also one; if f,-; is a Hurwitz polynomial and if f, has positive co- 
efficients, then /,, is a Hurwitz polynomial too. 
Another algorithm. Assume that f,(z) and f,(—2) do not have common zeros. Then 
two polynomials r(x) and ¢(x) with real coefficients and with no higher degree than 


n — 1 exist, satisfying 
fn(x)r(x) + fr(—x)t(x) = 2. (8) 
From this it follows that 


f(x) fa-i(—2) + f.(—2)fa.-i(z) = 2 with 2f,-.(x) = r(—2z) + Ua), (9) 
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the degree of f,_, being at most n — 1. Any other polynomial p(x) of any degree, which 
satisfies (9) instead of f,_, can be written as 


p(x) = fa-i(X) + 8(x)f.(2), (10) 


with a suitable odd polynomial s(7) = —s(—-2). Hence f,-, is the only polynomial 
satisfying (9) with no higher degree than n — 1. From (9) it follows that 


ga) = 1 (11) 


for the even component of the product f(x) = f,(x)f,-,:(—2); g(x) is definite on J’ and 
on J”, it is even definite on J. The product f(#) cannot vanish on J. Also its simplifica- 
tion F(x) cannot vanish on J, since the degree of F(x) is odd and F(0) = 1. Therefore 
f(x) and F(x) have the same number of zeros for Re x > 0 and for Re x < 0 according 
to Lemma 2. Lemma 3 is applicable to f, and f,_, . Thus, if f, is a Hurwitz polynomial, 
fn-1 is also one with proper degree n — 1. If f,-, with proper degree n — 1 is a Hurwitz 
polynomial and if f, has positive coefficients, f, is a Hurwitz polynomial too, and this 
is a consequence of (9). 

2. Details concerning the second algorithm of reduction. The second algorithm will 
be useful for the development of the formula announced. Some necessary details will 
therefore be developed. We assume f,(2) = a,x" + --- + a, to be a Hurwitz polynomial 
of proper degree n with real coefficients. As already stated, the polynomial /,,_,(2) 
defined by (9) is also a Hurwitz polynomial with real coefficients and with proper degree 
n — 1. The method leading from f,, to f,-, can now be applied to f,_, and so on. Thus 
we obtain a sequence of Hurwitz polynomials 

| Se Ae Le ek (12) 
with f, as a constant; f, has the proper degree k and real coefficients; any two adjacent 
polynomials f, , f,-, satisfy 

f(x) fe-1(—2) + fi(—2)f,-1(@) = 2. (13) 
It means no loss of generality to assume 
f,(0) = a, = 1; (14) 
(13) and (14) then lead to 
f,(0) = 1 for eB: Gos oe oe: i. (15) 
This in turn causes positive coefficients for all polynomials f, (see Sec. 1, Note). We 
increase all subscripts in (13) by 1 and subtract the new equation from (13); hence 
p(x) f.(—x2) + p(—2)f.(x) = O with p(x) = fisi(x) — fi-s(@); f(x) and f,(—x) have 
no common zeros. Therefore, 


fna(X) — fe-s(x) = Ces Bf, (2) for k=1,2,--+ ,n-I] (16) 
with a suitable constant 
Cpu oD. (16’) 
In addition to (16), we write 
fi(z) = 1 + e,2; ¢, > 0. (16”’) 
Regarding the positive constants c, , C2, --* , ¢, aS given, we can solve the system (16) 


with regard to f., --: , f, . We find: 








1952] INTEGRAL IN THE THEORY OF LINEAR SERVOMECHANISMS 299 


l+ce2z2 —!] 0 GO -- @ 6 
l cr —!l 0 0 0) 
0 1 e« —! 0 O | 
f(z) = ; k= 2,3 (17) 
| 
0 0 0 O- © 1 ex | 
This is a representation of all Hurwitz polynomials of proper degree k with f,(0) = 1. 


Vice versa all determinants (17) with coefficients c; > 0 give Hurwitz polynomials. 
Another representation may be given by means of the determinants 


C2 —] 0 0 O 
lL cx —1 0 0 | 
0 1 ena 0 0| 
| = F(z, ,¢, Cx) (18) 
0 0 O - »« | ef 
We can write then 
f(x) = Flac, , Ce, *°* 5 Ce) + PUeji ee Cay *** 5 Cady (17’) 


the right-hand-side showing the even and the odd component of f, . The functions (18) 
have imaginary zeros in 2 or real zeros in ix, which can easily be recognized as the eigen- 
values of a Hermitian matrix. The result about the zeros of the components of a Hurwitz 
polynomial with real coefficients is well known and has been found by E. J. Routh. So 
far this represents a minor application of (18). 

We are now going to develop another formula for f, where only the coefficients a; of 
the given polynomial f, appear. For this purpose we introduce the column-vector 


a by. a. = 0 for k>n and for k <0; 
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with 7 components, the matrices 
H Gis sa saeco cake ee PB nee me 
the so-called Hurwitz determinants 
Dy, = sign a, = 1, D, = , for i or oe (19) 


and the column-vecto1 


of the coefficients of the polynomial 
h Ca b, yi 4+ er h at +. 1. 
We then consider the polynomials 
rs 1)"*¢,(—az) fil) = Wait); kk = 0,1, +++, 0 20) 


with the two significant special cases 
From (1 and (20) follows for n — / a 


and we derive from (16’), (20’) and (21) that w,_,(2) has the proper degree n — k. 
1 3 


2 
This means: the product f,(2)f,(—z) does not contain the powers 2"*"”’, x 


a This is expressed by 

9.8. = —Qj.15. (22) 
There is only one polynomial f,_, of degree n — 1 according to (9). Hence there is only 
one solution 6,_, of (22) for k = n — 1, and this leads to D,_, ¥ 0. Let D,., ¥ 0; con- 


sequently the matrix 1 IS of rank k + 1, while the matrix (, >» Rena) consisting 
of all rows but the last of ,., is of rank k. This very matrix appears in (22), so only 


exists, and therefore D, ~ 0. Hence 


one solution of (22) for f, 
D, + 0; cm 1,2, ++" n= 1. 23) 
All systems (22) have only one solution 8, , and this belongs to 


i = * .— l 
. a Dr-1 vceot+ + ] 24 
f (2) = = Ay ; t= =e il, 24) 
- dD dD, 
a Qs he Q,-1 


The proof is clear. The coefficients of f, are positive. We have D, = a, > 0, apD,-,D;,' > 
0, D, = a,D,-, and thus, 


D.. > © for i= 





1952 INTEGRAL IN THE THEORY OF LINEAR SERVOMECHANISMS 211 


The coefficients c; in (16) are the quotients of the highest-power-terms of f; and f;-, . 
Therefore 


( a. q= Fe.’ he for k= 2,3, --- ,a. (26) 


The inequalities (25) form the well known Hurwitz criterion of stability. 


3. The formula for Y. Let 


m 


Plu, v) = DY axyu'v' (27) 


i,k=0 


be a polynomial of two variables u and v. Let y(t) and z(é) be two functions with con- 
tinuous derivatives p'y, p'z up to the order 7, k = m + 1. We then set 


m 


P*(y,2) = D> axyp'y-p*z. (28) 
We introduce Q(u, v) = (u + v)P(u, v). Obviously, 
| Q*(y, z) dt = | pP*(y, z) dt = P*(y, 2) (29) 
We consider the special polynomials 
Q,(u, f.(Wfi-1%) + fc oOfinm — 2; kom 1,2, +++ m. (30) 
From (16) it follows that 
Q,(u, v) = (u + vj Ses fe-s@&) + Qi-i(, v). (31) 
Therefore, 
Q,(u, v) = (u + v) DY caf: (2) f,-,(). (32) 


k=1 


We apply (29) to (32) with y and z as solutions of (2), i.e., f.(p)y = 0 and f,(p)z = 0. 


Hence 
2/ y(dzt) dt = — [ Q*(y, z) dt = : ot (CW fh i | (33) 
Ja k=l b 
with f*,(y) = fe-.(p)y and f#,(z) = fe-s(p)z. Setting y = z and a = 0, b = we find 
the announced formula 
2Y =2 [ y(t) dt = S olf? i(y)o)”. (34) 
“0 k=l 


We express c, and f,_, according to (26) and (24). We obtain 


ne dk —_ on ee = 1) “do 
2Y = aa;'g + > De Dit (35) 


rel 


Qox-1,k Qox-2,h ta Mi-1,h 


with the initial values g, as explained by (4). In this formula, squared linear forms of 
the g, appear together with the coefficients a; of the given equation. Formula (35) has 
already a form which makes it independent of the restriction a, = 1. It holds quite 


generally. 
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no 
_ 
bo 


In the special case go = 1, @ = G2 = *** = Qn-1 = O we find 
SY = 3c. (36) 
This sum can be easily computed from (16). Addition of all formulae (16) gives 
= n—1 
| Ps ae 2= -2+fothfhrt 2. Cfior — fi-) = 2 De, if 


or > 1c; = coefficient of x in (f, + fr-1). 


Therefore 


oY = a,_,a,' + D;':| a eel le Se See (36’) 
This formula too is not restricted to a, = 1. 
4. Two applications. 1) We set a) = a, = 1, which is no essential restriction. All 


other coefficients of f,, may be variable in order to minimize Y according to (36). This 
means, that the sum of all coefficients c; is to be minimized under the restriction ¢,c. --- 
c, = 1. An elementary calculation gives Min 2Y = n for c, = c¢. = -++ = ¢, = 1 with 


% = BL p-e 
9 a + 


Pe 


f(x) = 2" + (” 4 2” + 


/ 


Qe (” : 3° bo (” 2), by 


This formula can be proved by induction on n. 
(2) There are servomechanisms with an arbitrary input 6,({) and with a servo con- 


“< 


(37) 


trolled output 9(¢). The control depends on 

e(f) = @,(t) — 8; (t) 38) 
and shall make | € | as small as possible. According to the definitions given in [2], « can 
be called the regulated variable and @, the regulating flow. Let the servocontrol be of 


the proportional plus integral type, 1.e. 


AoA + 2,0 = —aeE — | a, edt (39) 
with constants a; > 0 for 7 = 0, 1, 2, 3. Combination of (38) and (39) gives 
do é€ + aye + ave + ase = —Q)0; — 4,0; . (40) 


Due to the integral in (39), e(¢) tends to zero with increasing ¢ if the right-hand-side 
vanishes identically and if 
D, = 4;@_ — Gaz > 0. (41) 
Now we consider the case 
6, 0 for i, <= @; 6g; = Cl for =e. (42) 


Then e(é) is a solution of the equation (40) made homogeneous. If the servomechanism 
is to start from rest at ¢ = 0, the initial values are 


€ 0) = @ = = 0; (0) = qq, = —C; €(0) = 72 = 0. (43) 
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Application of (35) leads to 


: a ay + a. a; 
2Y = CC =———.. (44) 
a,D. 
It is obvious that Y becomes smaller with increasing a, . Therefore a, should be made 
as large as possible. For practical reasons (saturation and overcontrol of amplifiers or 
the like) an upper bound for a, is given. With this in mind we minimize 2Y for a given 


a2 by variation of a, . Setting b; = a;/a, we find 
‘ ; , 2b; b, 2b,(b; le 
Maer - tS Seto (45) 
2b, b; 
bs = bi{(bi + i.” — by}. (46) 


This gives the best design with respect to the important case (42). 
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VARIATIONAL PRINCIPLES OF EQUILIBRIUM OF AN ELASTO-PLASTIC 
BODY* 


BY 
Y. YAMAMOTO 
University of Tokyo 


In this paper the author attempts to establish the general forms of the variational 
principles of equilibrium of an elasto-plastic body, and to make clear the relations among 
the principles presented previously [1]-[6]. 

§1 Fundamental principle. Consider an elasto-plastic body’ which is stressed by the 
surface traction F’ prescribed at each point on the portion Sy of the surface S of the 
body, the surface displacement v; prescribed at each point on the remaining portion 
S, of S, and the body force A‘ prescribed throughout the interior V of the body. The 
stress and strain distributions in the body are assumed to be given by o°’ and e;; re- 
spectively. Then we may select some incremental stress-strain law* to hold at each 
point in the body. We assume that this law is such that any small possible change of 


the stress-strain state satisfies the condition 
| 60° dée,; > 0°, (1) 


which is a generalized form of the so-called uniqueness condition [7]. 

Now, suppose that the increment AF’ of the surface traction F' on S, , the incre- 
ment Av, of the surface displacement v; on S, and the increment AK‘ of the body force 
K* throughout V are added gradually to the body. Then the resulting distribution of 
stress, strain and incremental displacement becomes o°’ + Ao’, €;; + Ae;; and Au, 
respectively. In the following lines we shall establish the variational principles which 
determine the resulting stress-strain state. 

Since the resulting displacement and strain are usually small, we apply the infini- 
tesimal deformation theory to our problem. An artificial displacement Au* , which is 
continuous and has piecewise continuous first and second derivatives, and which takes 
the prescribed value Av; on S, , is called an admissible displacement. Corresponding 
to it we may determine the strain Ae* by the equation 


4 
Ae¥, = 3(Au* ; + Au*,) ©. (2) 
Then the fundamental principle may be stated as follows: The following expression is 
non-negative for any artificial admissible displacement process Au*(f) (4 S ¢ S t,) 

*Received Oct. 29, 1951. 

‘We apply the Green’s theorem in this body, so that its surface must have a suitable regularity. 

2We need not assume one to one correspondence of the increments of stress and strain. 

3As usual the summation convention is used in this paper. 

‘Suffix ‘‘, j”’ is the sign of differentiation by the ordinate 2’. 
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such that Au*(f,) = Au; , the actual displacement in the equilibrium state, and Au*(t,) = 
Au* , an admissible displacement; 


* 


pae* 
il dV |  (o° + Ao**'(t)) dAe*(t) — | 


n € 


| dV p(K‘ + AK*‘)(Au* — Au,) 


me I dS(F' + AF")(Au* — Au,) °, 


ve 


S 5 


where p is the density of the material, the strain Ae*(‘) is derived by Eq. (2) from 
Au*(t), and the stress Ao*'’(t) is determined by Ae*(¢) and the stress-strain law. 
The expression (3) is equivalent to 


ap ppp 


IT dV(c"’ + Ao bAc*, — I dV p(K' + AK") 6Au* 


e 


(4) 


_— [| dS(F° + AF’) 6Au*® + I dV [ 6Aa*'’(t) dbAe*(t), 


vd 


where 6Au* = Au* — Au, , etc. By the principle of virtual work the first three terms 
in this expression vanish. Furthermore, the last term is always non-negative by Eq. (1). 
Thus, our proposition is proved. 

§2, Minimum principles. Unfortunately, the fundamental principle is not generally 
useful, because the value of the integral 

| ao + Ao*'’(t)) de *(t) (5) 

usually depends upon the stress-strain process. Hereafter we will assume that the value 
of this integral is determined only by the final values of Ac*'’ and Ae* at each point 
in v, and that any stress-strain state is attainable from an arbitrary state. In such cases 
the following principle is clear. 

Principle I: The expression 


ppp . 


U + AU* = || dV [ (o'’ + Ao*''(t)) dAe*(t) 


pee a 


_ II dV p(k’ + AK‘) Aut — | dS (F° + AF’) Au* 


Spr 
takes its minimum value when the admissible displacement coincides with the actual one. 
Certainly, the first variation becomes, 


ap 


6(U + AU*) = [{/ dV(a"’ + Ac*"’) 6Ae¥, — [| dV p(K* + AK‘) 6Au* 
- [| dS(F" + AF") 6Au* . 


‘Integration {> *‘* dAe¢*.(t) is taken along the process of the strain Ae* (lt). 
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From this we obtain the well known equilibrium conditions of stress as follows: 
(o'’ + Ao*"’), ; + p(K* + AK‘) = 0, in JV, (6) 
except on the discontinuity surface S, of stress in V, 

[(o,’ + Aot"’) — (o’ + Ao*'’)|]m; = 0, on S, , (7) 
where m, is the unit normal vector on S, , and o'} and o'! are the values of o*’ on the 
both sides of S, , ete., and 

(o'’ + Ao*'’)n; = F' + AF’, on S,, (8) 
where 7; is the unit vector in the direction of the external normal of surface. The second 
variation becomes 


6(U + AU*) = Ii dV 5Ac*''(t) ddAe*(t). 


0 
Since by Eq. (1) this value is non-negative, U + AU* takes its minimum value when 
the admissible displacement coincides with the actual one. 
By (2) we may easily obtain the equality 


i dV o'' Ac* — [{/ dV pK‘ Aut — [| dS F' Aut 


say. (9) 
= I| dS o*’n; Av; = const. = U, 
Then the following principle may be obtained by Principle I and Eq. (9). 
Principle I’: The expression 


AU* [// dV _ Ao*'’(t) dAe*(t) — [// dV p AK’ Aut — i/ dS AF’ Au; 


we 


Sr 

takes its minimum value when the admissible state coincides with the actual one. 
§3 Maximum principles. Since the incremental strain is given by Eq. (2), the in- 
volute transformation of Principles I and I’ may be easily obtained by the general 


method [8] [9]. 
In the first place, Principle I’ is equivalent to the following variational problem 


without any additional conditions, 


AU’ Il dV | Ao**"'(t) dAe**(t) 


_ [| dV p AK‘ Au*t* — [| dS AF’ Au** 


(10) 
za if dV n| acts _ * (Aut + Au**) 


— i/ dS p'(Au** — Av;) = Min., 


Se 
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where Au** is an artificial displacement which is continuous and piecewise differentiable, 
Ae** is an artificial strain which is independent of Au** and piecewise continuous, 
Ao**"’ is the stress corresponding to Ae** , \'’ and yw’ are Lagrangian multipliers, and 
especially \*’ is piecewise continuous and differentiable. Then the natural conditions of 
Eq. (10) become 


\ 


Ac¥* = 4(Aut* + Au**), in V,( 
(11) 
Au?" =.f9,.co &8;., 
Ag =} na V, i'n; = gw, on S,., 
AP = hn: . On Sr , 
\'; + p AK’ = 0, in V except on the discontinuity surface S, of \"’, and (12) 


(A; — AZ)m; = 0, on S,. 


Since, if we add Eq. (11) to Eq. (10) as the additional conditions, Eq. (10) returns 
to Principle I’, we substitute \"’ and y' into Eq. (10) from condition (12), i.e., 


apap f »de** 


aU’| =- I} aV( do* Act* — | 


vd “0 


Aa****(%) anet*(i)) 


id 


a I dS Aa**" nv; 


. (13) 
“ -[f/ dV . Ae**(t) dAo***'(2) 
+ I/ dS Ao***' n,; Av, = AU**, say. 
Now, the value of the integral | 
[~ Aci*(t) dAo**""(t) = Ao®*"" Act® — [~ Ao***''(t) dAet*(t) (14) 


is also determined by the final values of Ao**'’ and Ae** only. Accordingly, we may 
regard Ao**'’ as the independent variable in the expression (13). When Ao**"’ is piece- 
wise differentiable and satisfies the equilibrium conditions (6), (7), and (8), we call 
such a stress Ac**'’ an admissible stress. Then, by the general theory of the involute 
transformation the following principle is valid. 

Principle II’: The expression 


Ag** oa 


Ae*¥*(t) dAo**'’(t) + I| dS Ac**"' n; Av; 


AU** = — HH wad | ‘ 


Jde “0 


takes its maximum value when the admissible stress coincides with the actual stress. 


®Because the conditions (12) correspond to the equilibrium condition. 
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Certainly, we may confirm this principle as follows: Since, putting Ao**'’ — Ac’ = 
5Ao***’, it is clear that, 5Ac**'’n; = 5F' = 0, on S;, 


dAc**"’ = —6(pK') = 0, in V, 
except on the discontinuity surface S, of 6Ac**'’, and 
(dAc**"’ — §Ac**'’)m; = 0, on S,, 


we obtain the equality 


ap 


0 = I dV Au*t* 5Ac**") = — | dV Aut* 5Ao**! + I dS Aut* n; 5Ao**'* (15) 


for any function Au** which is continuous and has piecewise continuous first derivatives. 
Then the first variation of AU** becomes, 
sAU** = — I} dV Aet* 6Ac**!! + [| dS bAo**'' n,; Av, 


ep 


ee Hl dV(Aet* — Aut*) 6Ao**! + I dS 6Ao***' n,(Av, — Au**). 


That is, as the natural condition of this principle, we obtain the fact that Ae** may 
be derived by Eq. (2) from suitable possible displacement Au**. Furthermore the second 
variation of AU’** becomes, 


“0 


sause = — fff av| [sete dao**"(g 


pac 


_ Ae, ;(t) dAo'’(t) — (Ao**"’ — Ao"’) ac, | 


ep rb Aat* 
= -|// av | 5Act*(t) déAo***'(1). 
Since, we may attain to the stress-strain state (Ao'’, Ae,;) from the state (Ac**'’, Ae#*) 
by suitable process, we obtain the equalities 


rbidact** nbAdet* 


| 5Ac#*(t) déAo**"'(t) = SAo***? 5AcH* — 5Ao**''(t) ddAc**(1) 


“0 
»—bAe** 
‘ | 8’ Ao**"'(t’) dé’ Aet*(t’), 
where 
5’ Ac#*(t’) = Ac**(—t’) — Act. 
Accordingly, by condition (1) the value of 6°AU** is always non-positive. This shows 
that AU** generally takes its maximum value when the admissible stress coincides with 
the actual stress. 
By Eq. (9), we may obtain the following principle. 
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Principle II: The expression 


ppp »Ac** 27 


U + AU** = — Al dV | Ae**(t) dAo**'’(t) + I dS(a"’ + Aa***’) n; Av; 


vv “0 © 


Ae 


takes its maximum value when the admissible stress coincides with the actual. 
From the general character of the involute transformation the following relations 
are obvious, 
Min. AU* = Max. AU**, and 
(16) 
Min. (U + AU*) = Max. (U + AU**) = U + Min. au«s 


Certa inly, 


Min. AU* = AU* 


app side ppp ap 


on || dV | Ao''(t) dAe;;(t) — I dV p AK' Au; — I dS AF" Au, 


vee Jde vd 


SP 


2f Pp 


= I{/ dV Aco" Ae.; — || dV p AK' Au; — I/ dS AF’ Au; 


vide 


Aa ‘(t) ad«.,(0) 


= AU** = Max. AU**, 


+e Agi 
a 4 3 


§4 When are our principles valid? As already described our variational principles 
are valid when and only when the values of the integrals (5) and (14) are independent 
of the stress-strain process. Accordingly, when and only when the value of the integral 


a de; ; (17) 
] 


is determined by the final values of stress and strain only, our principles are valid. 
There are three reasons why the value of the integral (17) depends upon the stress- 
strain process: The first is the dependence of the stress-strain law upon the plastic 
history. Accordingly, we must assume that the stress-strain law is independent of the 
instantaneous plastic history. The second is the irreversibility of the stress-strain relation 
in case of the plastic deformation. That is, the stress-strain relations for loading and 
unloading states are individual ones. Then, our principles are valid for the cases where 
we may expect that the actual incremental process does not contain any unloading 
process at any yield surfaces except the initial one. Because in such cases we need not 
adopt the artificial admissible processes which contain such unloading processes. The 
third is the lack of the complete integrability of the integral (17) in stress space. We 
may easily obtain the general form of the stress-strain law which makes the integral 
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(17) completely integrable if the unloading processes do not occur, and it becomes as 
follows: 


’ kl Ta of pa of i 
des; = Eis do™ + F(f) +; (af + | df | 20% ——) , (18) 
da"’ 0a"" 
where F,;,,; = E,.;; is the elastic modulus, f(c) is the so-called loading function which 
is a function of c'’ only and the domain f < const. = c is the instantaneous elastic 


domain, and Ff) is a positive function of f(¢) only [10]. In this paper we assume such 
a stress-strain law.’ 

As above described, even if the stress strain law is given by Eq. (18), we can not 
assert the validity of our principles. But, if we can foresee that the stress-strain process, 
by which the resulting stress-strain state is attained, does not contain any unloading 
processes at any yield surfaces except the initial one f(o) = co , where c) is determined 
by the plastic history just before the incremental deformation, we may form our prin- 
ciples for the admissible states which are attained by such processes. If we assume the 
following reversible stress-strain law which satisfies the condition (1), we may form 
such principles, 


dAe = K..,, dAc™’, wherever f<a,orf =e& but df < 0, 
’ l ver of af l - p af rf 
Oa’ Oa 00°" 


wherever f > ¢ , or f = ¢ but df > 0, 
where f = f(o"’ 4. Ao’). 


The solution state of such principles satisfies the equilibrium and compatibility condi- 
tions, and by our assumption it is attainable by a suitable process whose stress-strain 
relation follows the Eq. (18), and which does not contain any unloading processes except 
those on the yield surface f = c) . That is, in this case the actual resulting state is gov- 
erned by these principles. 

§5 Principles for the stress and strain rates. Usually the mechanical quantities, such 
as strain, vary with finite time rates during a gradual stressing. In such cases they may 
be regarded to relate linearly with time during small time interval r, i.e., 

Ae*(t) = e*(), ete, (<t<7 
except in the very small region whose volume vanishes with 7. In the limiting case where 
r tends to zero, we may regard that the process, by which the resulting state is attained, 
does not contain any unloading process except those on the initial yield surface. Then, 
as already described, our principles may be formed as principles for the rates. That is, 
they are easily obtained from Principles I’ and IT’ as follows: 


The expression 


Me <= iii dV : ot ex — Il dV pK int - I/ dS Fiivt 


= vde 


Sr 


7This is a generalized form of the Hodge-Prager stress-strain law [1]. 
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takes its minimum value when the admissible displacement rate i* coincides with the 
actual one. 

The expression 


- 1 -: ok Bere ee ee a 
rr, : 5 U** = — [|| di = oer et * + I dS o**''n,v; 


~ Je 


takes its maximum value when the admissible stress rate ¢** coincides with the actual 
one. 

The principle of Hodge and Prager is equivalent to our Principle ITZ [1]. 

We may write down our principles for specific incremental stress-strain laws which 
belong to our general form (18). For example in the following lines we shall establish 
the principles for the material whose stress-strain relation follows the Prandtl-Reuss 


law,° 
P = (2G,,) - -- Ls 
where 
Bc l 
, ee “ar Iz = 588 
\0 wherever J. < k’, or J, = k* but J. < 0, 


and bu 
(2k°)~'s"’e,;; > O, wherever J, = k and J, = 0. 


In this case the following relations are easily obtained, 


€ == (). o "es = 2 € = (2G) 's 3 t- us''s’’, 
o'é;; = 2G,(€:;¢:; — ws''e;;) = (2G,) ‘3's’ > 02 


Then our principles become as follows: 
The expression 
7 


ui: 50+ = |] / av Genes — wts'e*) — [// dV pK‘ut — [| dS Frit , 


Jvdd 


» F 
where 
\9, wherever J, < k’, or Jz = k* but s'’e* < 0, 
1 * = < 
( (2K?) "Bes wherever J. = k and s'’e* > 0 
and the expression 
, ry F 4 Y ra ij? ij j ye F ° 
IT}: 2 l <-> — [|| dl ( 1G) lll cima + || dS o** nv; 


take on minimum and maximum values respectively when the admissible rates coincide 
with the actual ones. 


8The Prandtl-Reuss law and the Levy-Mises law belong to the Hodge-Prager’s general type as 
the extreme cases [1]. 
°This condition corresponds to condition (1). 
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The principles of Greenberg are equivalent to these principles [2] [3]. 

In the case where the elastic strain is very small compared with the plastic strain, 
we may neglect the elastic strain in the stress-strain law. In such cases our problem 
becomes somewhat complicated, because the stress state is not uniquely determined by 
a strain state. But the essential properties of our problem are unaltered by such cir- 
cumstances. For example, we derive the principles for the material whose stress-strain 
relation follows the Levy-Mises law” i.e., 


de;; = ds", 
where 
0, wherever J, < k’, or J; = k* but dJ, < 0 
ee 
2k*)~'s"’ detj , wherever J, = k’ and dJ, = 0. 


In this case the increments of displacement and strain increase linearly with time, but 
the stress change is independent of time, i.e., 


Ac® = e%(0), Ao**’ = const., etc., (0O<t< 7). 
Then the following relations are obtained, 


de;; = 8°’ [de,, de,.(2k’)*']'”’, 


v7 


| (o'? + Ao*''(t)) dAe*(t) = (0'' + Ao*"’) Ac* = (2k° Ac,® Ac*)'” 


~Aat* 


Ae¥*(t) dAo**''(t) = 0, 


| ba'’(t) dée;,(t) = 50°’ 5e;; = bs'’s'’[be,, 5€,,(2k")"]'” > 0. ° 


Then our principles are easily obtained by Principles I and II as follows: 
The expression 


I: U* = Ii dV(2k%,,é,,)'? — If dV pKint — I dS Fie 


nN 


Sr 


and the expression 


- 


‘.; U+* = I| dS(a"’ + Ao**'’)n 30; 


1 
take on minimum and maximum values respectively when the admissible displacement 
rate u* and the admissible stress Ao**'’ coincide with the actual ones. 

The principle of Markov and Hill and the principle of Sadowsky are equivalent to 
I. and JI, respectively [4]. 

§6 Conclusion. The fundamental principle and Principle I are valid even in the case 
of finite deformations, but the other principles have an essential restriction in condi- 
tion (2)."° 


10The Phillips’ second principle seems to be dubious in some respects [6]. Involute transformation 
for finite deformation becomes very complicated. 
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By the thermodynamical method we may prove the fundamental principle for stable 
equilibrium. Then the condition (1) seems to be essential in case of stable equilibrium.” 
Moreover, it is notable that usually the condition (1) is equivalent to the uniqueness 
condition. 

As described above our minimum and maximum principles are valid only in the 
restricted case. Accordingly, for the general cases we must apply the principle of virtual 
work directly, instead of the variational principle. 
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DIFFICULTIES WITH PRESENT SOLUTIONS OF THE 
HALLEN INTEGRAL EQUATION* 


BY 
S. H. DIKE 
Radiation Laboratory, The Johns Hopkins University 


I. Introduction. In a recent paper,’ rather serious discrepancies were shown to exist 
between values of broadside absorption gain and back-scattering cross section as found 
from experiment, and those predicted from Hallén’s first-order solution’ as modified by 
King and Middleton.* In this paper these discrepancies and certain additional short- 
comings with the present solutions to Hallén’s integral equation will be discussed. 

II. The first-order current distributions. It can be shown that the current distribution 
on a receiving dipole antenna is given by* 


I,(z) = Ig(z) — I,(2), (1) 


7 


where J,(z) is the current distribution due to the external field, EZ; , on the antenna 
with zero load (shorted), and J,(z) is the current distribution along the antenna when 
driven by a voltage, V; , equal to the voltage drop across the receiving antenna load. 
The distribution of Eq. (1) need be considered only in the cases where scattering be- 
havior is desired. All the other properties of a receiving antenna usually of interest, 
such as absorption gain, impedance, and effective length, are determined by the driven 
current distribution alone. 

The two different current distributions involved are given as follows: 

a. The transmitting dipole. The first-order solution of Hallén’s integral equation for 
the current distribution on a center-fed dipole of length 2h and radius a is given by 


I,(z) = j2nVof.(2)(tv.H2)', (2) 

where V, is the voltage at the terminals, ¢ = 1207 ohms, and 
f(z) = fil) + if’ = b, cos Bz — b, sin B | z | — C@) sin Bh + S() cos Bh, (3a) 
b, = [2y, + E(h)] sin Bh — S(h), (3b) 
b, = [2y, + E(h)] cos Bh — C(h), (3c) 
H, = Hi + jH?’ = [v. + E(h)] cos Bh — C(h). (4) 


*Received September 25, 1951. 

1S. H. Dike and D. D. King, The cylindrical dipole receiving antenna, Tech. Report No. 12, Radiation 
Laboratory, Johns Hopkins University, 1951. (Submitted to Proc. of 1.R.E. for publication). 

2. Hallén, Theoretical investigations into the transmitting and receiving properties of antennas, Nova 
Acta, Royal Soc. Sciences (Uppsala) 11, 1-44 (1938). 

8R. W. P. King and D. Middleton, The cylindrical antenna: current and impedance, Q. Appl. Math. 3, 
302-335 (1946). 

‘R. W. P. King, H. Mimno, and A. Wing, Transmission lines, antennas and wave guides, McGraw-Hill 
Book Co., New York, p. 163; 1945. 
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The functions C(z), S(z) and E(z) are defined in Ref. 3. The function y, is the ex- 


pansion parameter. It is determined in the King-Middleton method by considering a 


function defined by 


V,(z) = | g.(Z, 8)? 'e'”” ds (5) 
where 
r= [(z—s)? + a’]'” (6) 
and 
g.(z, s) = f(s)/f(z) (7) 
such that 
I,(z) = Io.f(z); I,(s) = Iof(s), (8) 


where J,, is the terminal current in the driven dipole. In the limit of vanishing dipole 
radius the driven dipole current distribution can be shown to be’ 


T,(z) = Ip, (sin Bh)~* sin B(h — | z]). (9) 
King and Middleton choose the function f(z) = sin B(h — | z |) giving 
g.(2z, 8) = sin B(h — | 8 })/sin B(h — |2]}), (10) 
and 
v,(2) = C(2) sin Bh — S(z) cos Bh (11) 
sin B(h — | z |) 


King and Middleton then argue that V,(z) is predominately real, and that a suitable 
expansion parameter may be found by setting 

V,(z) = | ¥,(2) | =v. , (12) 
where 2 is chosen so that W,(z)) is a good approximation to W,(z) over most of the 


antenna. Accordingly, they choose 


C(0) sin Bh — S(O) cos Bh | (sin Bh)™’," Bh <5, 


{| o(s —) sin ar — (1-3) cos h|, an > Z 


b. The unloaded receiving dipole. The first-order current distribution on a shorted 
electric vector of a plane-wave, far-zone field is 


(13) 


dipole antenna placed parallel to the 


T(z) = JArk fe z)(Biy,H,)” ee (14) 


where 
z) = 2y, cos Bz + E(z) cos Bh — C(z) + Cth) — cos Bh[2y, + E(A)], (15) 


H, = Hi + oot" = vy, cos Bh + E(h) cos Bh — C(A). (16) 


tromagnetic waves, D. Van Nostrand Co., New York, 1943, p. 142. 


5S. Schelkunoff, Elect 
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Fic. 1. The Function ¥,(z) for Bh = 0.4. 
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Fic. 3. The Function ¥,(z) for Bh = x 
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Fic. 5. The Function ¥,(z) for Bh = 2x 
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tbo 
aC) 
bo 


For the indefinitely thin shorted dipole, the current distribution is given by° 
Te(z) = Ipz(1 — cos Bh) ‘(cos Bz — cos Bh). (17) 


By the arguments of King and Middleton one is led to choose f(z) = cos Bz — cos Bh 
so that 


cos Bs — cos Bh 


qg,(z, 8) = : (18) 
9 cos Bz — cos Bh ) 
and V,(z) becomes 
C(z) — E(z) cos Bh 
¥,(2) == 19) 
cos Bz — cos Bh 
This function is plotted in Figs. 1 to 5 for various values of Bh and for 2 = 2 log 
2h/a = 10. By consideration of these figues and using the arguments of King and 
Middleton, one finds that a suitable choice of y, is 
C(0) — E(0) cos Bh | (1 — cos Bh)", Bh < a, 
y, = ] (20) 
e C(O) + E(0) |, Bh>r 


In the limit of very small Bh, y, = Q — 2andy, = 2 — 1. 
For the broadside case both the currents J,(z) and J,(z) are even functions. Equation 
(1) may now be expressed in terms of the foregoing as 


Ark; {, j2nZ.f (0) f (z)\ 
— fz2(2) -—— 


(z) = = ; (21) 
Bey.H, | fy, 
where Z, = £eb1/(Ze + Zr). (22) 
Z, is the antenna impedance defined by 
Z. = Vo/Iw .- (23) 
Z_ is the receiving antenna load impedance. For matched load 
Z,=iZ,\"/2R., Z,=Z. (24) 


III. Formulas for broadside gain, effective length and back-scattering cross section. 
Three equivalent expressions for absorption gain were derived in Ref. 1. They are 


gut Z. | | Bg.(h) (25) 
™ cViR, | H2 |? 
Bq,(h) , 
‘ oF 
G, Qy.T (26) 
' reid!" " 
CG, = 3 (27) 
ca 
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where 
T = Hy’ fi(0) — Hf"), (28) 
and where, for the broadside case,* 


Bg.(h) = 6 | f.(2) dz = (4y, — 20)(1 — cos Bh) 
™ (29) 


+ 2 Kin 28h + cos Bh(Eim 26h — Ein 46h — 4 Ein Bh) — j sin Bh Ein 46h. 


R, is the radiation resistance defined as the real part of Eq. (23). The symbol d denotes 
the effective length of the antenna. For the broadside case’ 


dui - [. 1.(e) de, ' (30) 
oo 
| d|/X = | Bg.(h) | (2x | f,(0) |)", (31) 
in 
d|/X = |Z, | | Bg.(h) | (fy. | He |)". (32) 
The relation for the back-seattering cross section, a, derived in Ref. 1 is 
o/d” = | Bge(h) — BZ.f2(0)Bg.(h) |? (rv? | Hy |)", (33) 
where 
B = j2n(ty,H.)' = In(HY + jH (ty. | He |)", (34) 
and 


Bgz(h) = B | fe(z) dz = sin Bh(4y, — 22 — 4 log 2 


j Eim 46h + 2 Ein 28h + Ein 48h) (35) 
+ cos Bh[2Bh(Q — 2y, + 2 log 2) — Bh Ein 48h — 26h Ein 26h 


— j Ein 46h — 2sin 28h — 2j(cos 28h — 1)]. 


IV. The required equality of y, and y, . Although the method of King and Middleton 
yields a different expansion parameter in the receiving case from that found by them 
for the driven dipole, it is necessary, in order to have a consistent theory for the re- 
ceiving dipole, that both expansion parameters be identical. This can be seen from a 
comparison of the two equivalent definitions for effective length. In addition to Eq. (30), 
effective length may be defined by 


d = I,pZ,/E; . (36) 


*The function Ein(z) is defined in the Appendix. 
1S. H. Dike, The effective length of antennas, Tech. Report No. 13, Radiation Laboratory, Johns 
Hopkins University, 1951. (Submitted to I.R.E.). 
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From (30) and (36) one obtains the interesting equality that 


a a 
E77, |, 10%. (37) 


Substituting in (37) from (2), (14), and (29), one obtains 


2f2(0) _ Bg.(h) (38) 


YH, y.H,° 
This relation is true if Y, = y, , because when this is so, H, = H, and 2f,(0) does indeed 
equal 69,(h) for all orders of solution. A second relation is possible between y, and y, 
which satisfies (38), but this relation is dependent upon the order of the solution. This 
means that the expansion parameter would be a function of the number of terms retained 
in the solutions for the currents. This is obviously undesirable. In any case y, must 
equal y, for the zero order solution. 

This fact constitutes the first difficulty encountered in the King-Middleton method. 
It is difficult to say which of the two parameters, y, or y, , is the better. It appears 
that neither is particularly good.’ 

V. The behavior of the theory for short dipoles. It was pointed out in Ref. 1 that the 
yalues of gain obtained from Eq. (26) do not reduce in the limit of decreasing Bh to the 
value 1.5. One would expect the value 1.5 as being the correct one for finite 2 because 
for very small 6h the current distribution must be essentially linear. This does not 
imply that either the radiation resistance or the effective length should reduce to those 
of the indefinitely thin short dipole. For very small 8h, Eq. (26) becomes 


3 (2y, — 2+ 3)° ; 
ee ee 2y.(Q — 2 — 2 log 5} oniaands 32) 


This result is independent of Bh but remains a function of the expansion parameter. 
Since for small Bh, ¥, = Q — 2, Eq. (39) becomes 


3 (Q — 1)° 
G, = | ——— “+ h<1. (40 
"= sl@—DAN2-2+4logny "MS 
This is not 1.5 except for Q ~o. It is not evident that the difficulty would be 
removed by retaining more terms of the series solution. Various values of short dipole 
gain can be obtained from (39) by using the expansion parameters of previous authors. 
This is shown for 2 = 10 in Table I, where the value of the expansion parameter for 


small Bh is given. 


TABLE I 
Author yp G, Percent Error 
Hallén 2 1.5114 0.76 
Gray® 2 — 2+ log 4 1.4835 1.10 
King-Middleton 2-2 1.4098 6.01 


If the parameter y, is used, where y, = 2 — 1 for Bh < 1, the value of G, at 2 = 10 is 
1.464 which is an improvement over the use of y;, . 


8M. C. Gray, A modification of Hallén’s solution of the antenna problem, J. App). Phys. 15, 61-65 
(1944). 
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The fact that the bracket of (39) should be unity in the first-order theory gives 
VY — (22 —8 + 4 log 2) + & —62+9 =0, 
or 
Y = 2 — 2.6138 + 0.879 (Q — 2.807)'”. (41) 


Of the two values allowed, the larger is probably the one that should be chosen. 

A similar situation results in the value of o/\” for very short matched-loaded dipoles. 
The theory should reduce to the value 9/167. In the limit of decreasing Bh, Eq. (33) 
for the matched-loaded case becomes 


o 9 | ¥:. — 2 + 3)(2y, — 2+ DY (42) 


\? 16% 1 — 29- 





V(3¥, — 22+ 44 4 log 2) 
If the requirement is made that y, = y, then 


i. E __(2y — 2+3) _ | 
” 16r LBV” — 2y¥(Q — 2 — 2 log 2) ]° 
The bracket of (43) is identical to that of (39) leading to the same requirement on 
given by (41). It appears then that an additional requirement should be imposed on 
the expansion parameter that has not previously been considered. 
It is also of interest to consider the value of impedance in the limit of very small Bh. 
The impedance is given by 


(43) 





Z. = Vo/In. = — itv. H[2rf,(0)]"’, (44) 
which for Bh « 1, becomes 
7 Rp 2x — —WI3dY. + 2i(6h)*) 
2. = Re + jXe = sahoy, + x) + iGW)” (45) 
where 
x=2+2log2 —Q. (46) 
Separating real and imaginary parts of (45): 
R, = 20y,(8h)*?(3p, + 2x)(2p, + x)~’, (47) 
and 
oa (48) 


~ BR(2y, + 2)’ 
Since y, = 2 — 2 for Bh K 1, 








_ ___ 20(8h)*___ 
mes + (2 log 2)*/y’ (49) 
where 
y = QD + (log 2 — 1)(log 2 — 1 + 49), (50) 
and 
_  —60 (Q — 2)? | 
x. = Bh Ee ‘ (51) 
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For a value of ¥ which satisfies (41), 2, becomes 
R, = 20(Bh)?(2Qy — 2 + 3)?(Qyp — 2+ 2 + 2 log 2)~”. (52) 


Using w from (41) for 2 = 10, R, = 20(h)*(0.941). The corresponding value from the 
King-Middleton relation (49) is R, = 20(6h)*(0.978). There is a difference of about 
four percent between these two values. 

The corresponding values of X, for Q = 10 are X, = — (7.375)60/6h for y satisfying 
(41), and X, = —(6.819)60/h for the King-Middleton expression (51). The difference 
here is about eight percent. 

It is of interest to compare these results with the reactance computed from the static 
capacitance between two cylinders placed end to end in air, and separated by a distance 
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which is negligible compared with their individual lengths. This reactance is given by’ 


; —¢(Q — log 12) —60(Q2 — 2.485 7 
Ag == Pe NE gy — ee 9 (53) 
2nGh Bh 
For Q = 10, Eq. (53) gives X,, = —(7.515)60/8h. This result may be compared with 


those obtained from Eq. (48) listed in Table II. 


TABLE II 


—BhX. . 
y ar for2 = 10, Bh «<1 
60 
2 (Hallén 7.470 
from Eq. (41 a 7.375 
Q — 2 + log 4 (Gray) 7.246 
Y, ees 7.114 
Q — 2 (King-Middleton) 6.819 


VI. The use of the King-Middleton values of impedance. It was thought that im- 
proved values of the various antenna properties might be obtained by using the second- 
order impedance values of King and Middleton in the expressions for gain, back-scattering 
cross section, and effective length. Figure 6 shows the absorption gain, G, , computed 
from the straight first-order theory, and also as computed from Eq. (25) using second- 
order values of | Z, |"/R, , and from Eq. (27) using first-order effective length and 
second-order R, . As can be seen from this figure, very different results are obtained 
when the King-Middleton values are used. The use of Eq. (27) with second-order R, 
would seem to indicate that the King-Middleton values for the radiation resistance are 
too large in the region of Bh = 2, and are too small in the region of Bh = 3.5. The use 
of Eq. (25) with second-order | Z, |*/R, yields results which are certainly contrary to 


fact near 8h = 1.5 and Bh = 4.5. The “bulge” in the first-order theory near Bh = 2 is 
contrary to experimental data and disappears if y, is used instead of y,.’ 

Figure 7 is a comparison of the first-order theory for effective length given by Eq. (31), 
the relation of Eq. (32) using second-order Z, , and Eq. (27) using second-order R, and 
first-order G,. The use of (32) with second-order Z, gives results which are unreasonable 
near Bh = 1.5 and Bh = 4.5. 

Figure 8 shows the back-scattering cross section for matched load according to the 
first-order theory and also from Eq. (33) where second-order King-Middleton values 


are used for Z, . The latter curve behaves strangely near Bh = 1.5. Neither curve 
represents experiment, particularly in the region above Bh = 4.’ 


These three figures show that it is not permissible to use second-order King-Middleton 
values of impedance in first-order formulas. This may be due to the unknown behavior 
of the series solution as regards convergence,’” or it may be that the second-order im- 
pedance values of King and Middleton are not good. This latter case would imply 
that the expansion parameter y, can be better chosen. 


VII. The problem of choosing the expansion parameter. The integral equation of 


*R. W. P. King and C. Harrison, Jr., The impedance of short, long and capacity loaded antennas with a 





critical discussion of the antenna problem, J. Appl. Phys., 17, 170 (1944). 
108. Schelkunoff, Concerning Hallén’s integral equation for cylindrical antennas, Proc. I.R.E., 33, 872 


(1945). 
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Fig. 8. Back-scattering cross-section for Q = 10, Z, = Z,* 


a First-order theory 
——— Equation (33) using second-order | za |?Ra 


Hallén is known to be a sufficiently accurate formulation of the problem. It has been ex- 
amined by many workers,"’''*’'* and has been shown to contain approximations only to 
the order (a/h)’."° Hallén proposed an iterative process for solving this equation and ob- 
tained a series solution.” Modified solutions have been proposed by Miss Gray* and by 
King and Middleton.* These have consisted essentially of modifying the expansion 


1}). Middleton and R. W. P. King, The thin cylindrical antenna: a comparison of theories, J. Appl. 
Phys. 17, 273-284 (1946). 

2],, Brillouin, The antenna problem, Q. Appl. Math. 1, 201-214 (1943). 

188, Schelkunoff, Antenna theory and experiment, J. Appl. Phys. 15, 54-60 (1944). 
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Fic. 9. Back-s tering cross-section at Bh = 27 for 2 = 10, Z, = 0, as a function of 


expansion parameter 


parameter used. Hallén™ claims that such modifications have no mathematical founda- 
tion, and this may be so. However, the method of King and Middleton does appear 
reasonable from the standpoint of physical reasoning. The choice of a trial function 
which is known to be representative of the current distribution in the limit of vanishing 
dipole radius seems to have some merit, although the arguments as to why this should 
yield a better solution have been attacked by Hallén."” Nevertheless, if it is thought 

ME, Hallén, Admitta jra for antenn ind the relation bet 2 antenna theories, Tech. Report 
No. 46, Cruft Laboratory, Harvard Univ., 1948. | 

16E), Hallén, 7'ravel? and unsymmetrically fed antennas, Tech. Report No. 49, Cruft Labora- 


tory, Harvard Univ., 1948 


1952] DIFFICULTIES WITH SOLUTIONS OF HALLEN INTEGRAL EQUATION 241 


that such a choice of a trial function is justified, it appears from the foregoing discussion 
that consideration must be given to both Eqs. (11) and (19) in making a final choice 
of y. Note that these two expressions are identical at all resonant lengths (8h = 7/2, 
37/2, 54/2, ---). Also since ¥(z) is not predominantly real for the longer lengths, the 
final choice of Y may best be a complex value. 

Back-scattering cross section appears to be a property of the dipole antenna which 
is particularly sensitive to prediction by theory. As an illustration of its sensitivity to 
the choice of expansion parameter, lig. 9 shows the back-scattering cross section at 
Bh = 27, 2 = 10, for the shorted dipole as a function of the expansion parameter. A 
factor greater than two exists between the result using the King-Middleton parameter 
of about 6 at this length, and the Hallén parameter of 10 in the first-order theory. 

VIII. Conclusion. In view of the fact that the series solution of the integral equation 
has been studied, criticized, and modified by many authors since Hallén’s first paper in 
1938, and since a theory which can be practically computed does not seem to exist 
which adequately predicts the complete behavior of a simple dipole antenna, it appears 
perhaps that a new attack on the problem is justified. 

It is significant that the results of Van Vleck, et al,’® for the back-scattering cross 
section of a shorted dipole agree more closely with experiment than the first-order 
solutions of Hallén, King and Middleton, or Miss Gray. Such a comparison is made 
in Fig. 16 of Ref. 1. Hallén’s recent solution’’ for the driven dipole may be an improved 
one from the standpoint of gain. Some attempts have been made to solve the integral 
equation by variational methods.'*''* Storer’s solution fails for Bh greater than 32/2. 
Tai removed this difficulty but his first-order values of R, at the first resonant length 
are still higher than those of King and Middleton. It may be worthwhile to follow 
up a suggestion made by Brillouin’ that the known function and the kernel of the in- 
tegral equation be expanded in Fourier series with known coefficients, and that the 
unknown function for the current be likewise expanded with unknown coefficients. 
Term-by-term integration would then lead to a set of simultaneous equations for de- 
termination of the coefficients. No published results of such an approach have come to 


the author’s attention. 
APPENDIX 
Kin (x) = Cm (2) + 781 (a) 
az vat~ 
“i 2 cos ft 
Cin (x) = | dt 
Jo l 
re * gin ¢ 
sl (zr | di 
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ON BROWNIAN MOTION, BOLTZMANN’S EQUATION, 
AND THE FOKKER-PLANCK EQUATION* 


BY 
JULIAN KEILSON 
AND 
JAMES E, STORER 
Cruft Laboratory, Harvard University 


Abstract. In order to describe Brownian motion rigorously, Boltzmann’s integral 
equation must be used. The Fokker-Planck type of equation is only an approximation 
to the Boltzmann equation and its domain of validity is worth examining. 

A treatment of the Brownian motion in velocity space of a particle with known 
initial velocity based on Boltzmann’s integral equation is given. The integral equation, 
which employs a suitable scattering kernel, is solved and its solution compared with 
that of the corresponding Fokker-Planck equation. It is seen that when M/m, the mass 
ratio of the particles involved, is sufficiently high and the dispersion of the velocity 
distribution sufficiently great, the Fokker-Planck equation is an excellent description. 
Even when the dispersion is small, the first and second moments of the Fokker-Planck 
solution are reliable. The higher moments, however, are then in considerable error—an 
error which becomes negligible as the dispersion increases. 

1. In the treatment of Brownian motion, it is customary to assume a Langevin 
equation and simple dynamical statistics of the individual collisions and then to deduce 
a Fokker-Planck equation describing the random motion of the heavy particle. The 
Fokker-Planck equation obtained is a second-order partial differential equation and the 
absence of higher-order differential terms is inferred directly from the above assumptions. 
As will be seen, the solution of this Fokker-Planck equation does not provide a com- 
pletely satisfactory physical description. Consequently, the assumptions underlying the 
equation cannot be correct [1, 2, 3, 4] and the extent of their approximate validity comes 
under question. 

That the solution of the Fokker-Planck equation is not a wholly satisfactory repre- 
sentation of Brownian motion may be seen in the following way. Consider a heavy 
particle known to have the velocity v) at 4 = 0. For all subsequent time, there is a finite 
probability that the particle will have undergone no collision. It must, therefore, be 
expected that the probability density w(v, t)** describing the stochastic motion in velocity 
space will always have a singular component of the form f(/)6(v — vo), where 6(v — vo) 
is the Dirac delta-function. If one were to try to describe the motion by the Fokker- 
Planck equation 


aw(v, l 2 | 
= “ we 5 DV ‘uly, t) + 1V-{vu(v, 0}, ) 
F 2 


*Received Nov. 3, 1951. The research reported in this document was made possible through support 
extended Cruft Laboratory, Harvard University, jointly by the Navy Department (Office of Naval 
Research), the Signal Corps of the U. S. Army, and the U. 8S. Air Force, under ONR Contract N5ori-76, 
1.1 

**The function w(v, 0) obeying (2) is often represented in the literature by P2(vo/v; t), the probability 
density for velocity v, ¢ seconds after there is a known velocity v. 
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subject to the initial condition 


wiv, 0) = &(v — v,), 2) 


no such singular component would be available in the solution. The immediate dis- 

appearance of an initial singularity is, indeed, characteristic of all diffusion equations of 

finite order. Only by means of an integral equation can such a singularity be maintained. 
All of this, of course, is in keeping with the fact that fundamental to the description 

of Brownian motion is Boltzmann’s equation, an integral equation of the type desired [4]. 

If A(v’, v) dv is the probability per unit time that a particle with velocity v’ will undergo 

a transition to a volume dv about v, the Boltzmann equation describing the motion is 

tate 4 ; , 


w(v’, t)A(v’, v) dv’ — wiv, ¢ A(v, v’) dv’. (3) 


ol 


ry 


This is simply an expression of the fact that the rate of change of the population of a 
cell in velocity space is the difference between the rate of departures from the cell and 
the rate of arrivals. 

From this Boltzmann equation a corresponding Fokker-Planck equation may be 
derived. If Eq. (3) is multiplied by an arbitrary, but suitably behaved function R(v), 
and integrated over v, 


e af 


Ow(v, t) , 
| R(v) = rv dy = iT] R(v)w(v’, t) A(v’, v) dv’ dv 
_- I R(v)w(v, t)A(v, v’) dv dv’ 
= I | 2, = =v | — rw’) t) A(v’, v) dv dv’ (4) 
_ // R(v)w(v, t)A(v, v’) dv dv’ 
= [| pa = “d “we Eee h v’, t)A(v’, v) dv’ dv 
Here v-—v’) -V’ is to be understood as the dot product of two n-th-rank tensors, 
Integrating by parts, one has 
/ Riv’ » ae dy’ = I Riv’ >= Ve" fiw’ — vw) A(v’, v) Ww’, t)} dv’ dv. (5) 


0 


Since R(v’) is an arbitrary function, the associated coefficients may be equated to yield 


Ow(v, ft . ] 
> = VV": fA, (v)wlv, 1}, 6 
dt — ni 
where A,,(v) is the 
1,.(v | v-—v i(v, v’) dv’ 7) 
Equations (3 nd (0) are equivalent and provide an exact description of Brownian 


motion. 
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This treatment may be readily generalized to include Brownian motion in coordinate 
and velocity space. 

2. In those treatments of Brownian motion based on Langevin’s equation, moments 
higher than the second are found to vanish, and the Fokker-Planck equation (3) is 
obtained. As already noted, such an equation is certainly unsatisfactory when the dis- 
persion is small. It would, therefore, be desirable to try to treat the Boltzmann equation 
directly. Plainly an exact kernel A(v, v’) is unavailable and its use is almost certainly 
not feasible. However, it is possible to introduce a kernel which provides a reasonably 
accurate description of the microscopic scattering process and which is, at the same time, 
amenable to treatment. Such a kernel is of the form A(v, v’) = @(v’ — yv), where y 
is a dynamical damping parameter close in value to, but less than, one. Some justification 
for this form may be found along the following lines: 

Let B(v, v’) dv’ be the probability per unit time of a particle with initial velocity v 
making a transition to a volume element dv’ about v’, when all the particles with which 
the heavy particle collides are stationary. If the lighter particles have an equilibrium 
distribution w(v’’), then 


. 


AG, VW) = | Biv —v"’,v 


. 


, 


— v’’)w(v’’) dv’’. (8) 


Since the particle under observation is very much heavier than the particles with 
which it collides, B(v, v’) is a highly localized function of v’, centered roughly about 
yv where again y is very close to but less than unity. If B(v, v’) is assumed to have the 
form B(v’ — yv), then 


. 


A(v, v’) = Biv’ — v”’ — y(v — v’’)|w(v’’) dv’’ 


= | Biv’ — yw —(1 — y) v’’ |w(v’’) dv’, 
so that this will have the form of @(v’ — yv). 
Note that the form of @(v’ — yv) implies that the mean free time 7 of a heavy particle 
is independent of its velocity, since 


= = A(v, v’) dv’ = Q(v’ — yv) dv’ = | a(v’’) dv’’, (9) 
a constant. This behavior is proper to Brownian motion where the heavy particle moves 
so slowly compared to the lighter particles that the mean relative velocity of the heavy 
particle does not vary significantly. 
It would appear offhand that the functional form of @(v) could be chosen arbitrarily. 
However, this is not the case since @(v, v’) must satisfy the equilibrium condition: 


wv’) A(v’, v) = w(v)A(v, v’) (10) 


where w(v) is the equilibrium distribution of the heavy particles which the particle will 
ultimately assume. If it is also demanded that w(v) depend only on | v |, the two re- 
strictions imply that @(v) must have the form @, exp |—8v°} and that w(v) must 
have the corresponding form w) exp {—8(1 — y°)v’}, where @» and w) are constants 
see Appendix 1). That the Gaussian character of the equilibrium distribution follows 


from the form of @(v — yv’) is reassuring. 
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Thus, the Boltzmann equation to be solved is 


dw(v, t) j 2 ] ; 
— Qo | w(v’, t) exp {—B(v — yv’)"} dv’ — — w(y, ), (11) 
c J tT 


where 


| - , 
- = | Q(v’ — yv) dv’ = @, | exp {—B(v’ — yv)*} dv’ 


a 


3/2 


7 3 
= @, | exp {—v’”} dv’’ = @, () 
| exp t-av?} : 


Before discussing the solution of this equation, it is worth while to put down the 
corresponding Fokker-Planck equation. The first moment will be given by 


Ay =e | (v — v’) exp {—B(v’ — yv)’} dv’ 


= @ | l(yv — v’) + (1 — y)v] exp {—B(v’ — yv)*} dv’ 


(13) 
(l—y)v exp {—6(v’ — yv)*} dv’ 
a \?/2 { oe 
= a,(™] (1 — y)v = (7 “y. 
B/ T 
Similarly for the second moment, 
A; = Q | (v — v’)(v — v’) exp {—B(v’ — yv)°} dv’ 
= Q, [(v’ — yv)(v’ — yv) + (1 — y)’vv] exp {—B(v’ — yv)*} dv’ 
(14) 
Qo 3 si 
=— 7 Be + ( — wa) 
2 
1 — y)’vv 
—_ € + - —. a 
B Br 


where ¢ is the unit tensor. 
If the latter part of A, is ignored since (1 — y)* is small and if the higher moments 
(whose effect will be small for ¢ > 7) are ignored, Eq. (1) is regained where now 
] 1-—¥ 
D = — and ;=— ; (15) 
2Br T 
As is seen in Appendix 2, the solution of the Boltzmann equation (11) subject to 
condition (2) is given by 


\ > ( ’ 
, A r B \ t > 
walv, 0 -| v= Vor a 1 (4) (4 } ee es ee | exp {~! [> (16) 


. 


<i 
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where 
A on i-7° (17) 
_ Soa 1 pete 7 ° 
The solution of the Fokker-Planck equation (1) subject to condition (2) may be 
taken directly from Wang and Uhlenbeck [1] and is given by 


n “ —n(¥ — Vo exp Laity 

_— j= > Seine P | 
Weel, Y EF — exp a | exp | DQ — exp {—2nt}) (18) 
Note that the singularity 6(v — v,) is preserved in the solution of the integral equation 


but is not in the solution of the Fokker-Planck equation. For t > 7, however, the delta- 


function ceases to play a prominent role. 
From Eqs. (16), (17) one finds that the equilibrium distribution for the solution of 


the Boltzmann equation is given by 


oe _ awa 4 2,2 
w,(v) = lim w,(v, t) = “ exp {—@(1 — yu}. (19) 
For the Fokker-Planck equation, 
“ ” 3/2 . _ 1» 
Wrp(V) = (2) exp { D’ \ (20) 
Inserting the values of 7, D from (15), this becomes: 
insti = (20 — m6) exp {—28(1 — 7)”}. (21) 


Plainly if y is sufficiently close to one, then 
(l—y7’) = (1 —y(l + y) > 21 — ) 


and the two equilibrium distributions are identical. 
It is also of interest to compare the manner in which the average velocity and the 


variance vary in time. These quantities are defined by 
{v)(t) = [ vow, t) dv 


and 


. 


a(t) = | (v° — (v)*)w(v, t) dv. 
v)z(0), (v)rp(t), of(f), and opp(t) may be computed directly from their respective equa- 
tions of motion. Thus, if Eq. (1) is multiplied on both sides by v and integrated over v, 


then 


a \ 1 —, f 
AY = —7n(v), so that (v)pp(t) = Vo exp {—nt} = vo exp {t=}, (22) 


Similarly, multiplying by v° and integrating, it is found that 


“ = 3D — 2n\v’), (23) 
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which yields in turn 


do d 


di dt’ 


(v’) — (v)’) = 3D — 2no’; (24) 


and SO 


3D 3 2 \) 
Ovepl(t) = ( _ —— — . — av et. — Or 
rrp(t - I exp { —2nt}) 131 —(1 exp { ~ (I yt? }. (25) 


t 


The same procedure may be applied to the integral equation (11) to give 


f fi—y)\) 
(v)5(t) = Vo exp 4 ~( ah (26) 


and 


3 | ft =o") ez at . 
v'),(?) = 5 l— exp 4—( —4 )eh + Vv, exp \—( Y Je}. (27 


28(1 sole: jie 


Correspondingly, one finds that 


- 
én) = > | | — exp \— 


+ v,| exp me = U) — exp {-2(3 — 1), 


These same results could also have been obtained from the solutions (16) and (18), 
but the computations are more tedious. 

It is seen that (v f) and (Vrp)(t) are identical and that op(t) and Cer(t) are nearly 
identical. Indeed, if the smaller term in A. had not been ignored in obtaining the corre- 
sponding Fokker-Planck approximation, o;(¢) and of»(/) would have been precisely the 
same. For consider the Boltzmann equation in its differential form: 

ow 


> Vv -(A,w). 


ot ad 
The above procedure yields 


= / vV --(A,w) dv + [ vV~-(A.w) dv, 


since all integrals involving higher moments vanish when integration by parts is carried 
out. Moreover, from the choice of A(v, v’), the two integrals are simple functions of 
v) and the abov differential equation does determine v(t). The same procedure 
applied to the Fokke r-Planck equation can only vield the same result, because all con- 
tributing terms are present 

It is seen then that the validity of the Fokker-Planck approximation is excellent 
when ¥ is sufficiently close to one. For the ordinary domain of Brownian motion this 
will certainly be the case. For the elastic collision of hard spheres, for example, it is 


easily found that 


} m 


(sy) = — Vv 
3 M+m 
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where (év) is the mean change in velocity suffered by a particle of mass M and velocity 
in a single collision with particles of mass m. Then 


div) (ov os ~§ mm y 
dt +r 838 (M+m)r 


so that, from Eqs. (15) and (22) 


_d a oo | Mm 
9 “3 (M+ m)r 
| — y then is given by 4/3 m/(M + m) and for typical Brownian motion will be ex- 


tremely small. 

If it were possible to treat the exact kernel A(v, v’), one would still expect to find 
excellent agreement between the Fokker-Planck and Boltzmann solutions for ¢ > r. 
Even when ¢ ~ 7, the first and second moments of the Fokker-Planck equation should 
be reliable. But for  ~ 7, higher-order moments would be in considerable error. However, 
for { > 7, these errors will become entirely negligible. 
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Appendix 1 


The Restrictions on @(v) and w(v) Imposed by Equilibrium 


Denoting the rectangular components of v by v; , 7 = 1, 2, 3, and letting 
V0, , Vs , Vs) = In a), 


the equilibrium relation (10) can be written in the form 


(1-1) 
In w(v) + Wei — Yr , v2 — Yo , V3 — Ys). 
Taking the partial derivative of Iq. (1-1) with respect to v, , v/ , and noting that 
fa) . a 
——; Inw(v’) = =——x7 Ino(v) = 0, 
Ov, Ov Ov, OV 


it is seen that 


au # 


’ Pua atin — eee i es hats i Trt : 
YVs;(0U1 — Yi ,? yU2 , V3 — YU3) Yi YY 5 Ys YV2 , U3 YV3); (1-2) 
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where 
a” 
w= Ss = (V, , Vo , V3). 
v ] Ov, ov, v 1 2) 3 
Setting v} = yu; , v?’ = (1 — y’)v; in (1-2), it is further observed that 
v,;(v'’, v2’, vf’) = ¥,,(0, 0, 0) = —B.; , 


where 8;; is a constant. Hence W(v, , v2 , v3) must be of the form 


3 3 3 
Y(X; , 2,3) = — p 2 7. B,; vv; + z. av; + A, 
i 1 7=1 1 


where a; and A are constants. Inserting this result into (1 — 1), one finds that 
3 3 3 
Inw(v) = — >> 8,11 — Yow; + Dal + yu + A’. (1-3) 
?=1 7 1 1 


But, since the distribution of small particles is assumed to be isotropic, one has 


wv) = wo! Vv!) = w((v, + v3 + 05)'””). (1-4) 
The only possible way (1 — 3) can satisfy this condition is for 
B;; = Bé;;, a; = 0. (1-5) 
Hence 
Q(v) = @y exp {— Bu}, wv) = wo exp }—B(1 — yw}, 


where @, and w, are constants. 


Appendix 2 
Solution of the Boltzmann Equation 


It is desired to solve Eq. (11) subject to the condition (2). Two methods will be given. 
One procedure is to introduce the Fourier transform of w(v, ¢), 1.e., 


7 Ae...) = | exp {ik-v}w/(v, t) dv 
with 


T’,.(k, 0) = | exp {ik-v} 6(v — vo) dv = exp {7k-vo}. 


Taking the Fourier transform of Eq. (11) one obtains the following equation for 
7. (e:; f): 
6] 


: eas es : 
T (kt, t) = A(k)T. (yk, t) — — T.(k, 0), (2-1) 
ot T 
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where 


d 


_ (z)" p\—*} = = {- 
= Gol 5 ex] 46) ~ 7 xp 4a)" 


It is now convenient to introduce the Laplace transform 


A(k) = exp {7k-v}@(v) dv = Q@, [ exp {ik-v — Bu} dv 


L.(k, s) = | exp {—st}TW(k, 0) dt; 


“oO 
then 


| exp {—st} = T..(k, t) dt = —T,,(k, 0) + sL,(k, s) = —exp {ik-v,} + sL,(k, s). 


Thus, taking the Laplace transform of (2-1), one obtains the equation for L,,(k, s), 
k* l 
—exp {ik-v,} + sL,,(k, s) = — exp § —=/L,(yk, s) — — L.(k, s). 
T q 48 c 


This may be rearranged to give 


L.(k, s) = exp {ik-v,} + a | = exp {— En, (yk, s). (2-2) 


a are Ts+r 


The finite difference equation (2-2) may be solved by the following procedure: Replace 
k by vk This yields 
l : l l Ks 
L,. (yk, s) = + ; exp foyk-v)} + - — exp {—? Hh (y'k, s). (2-3) 
s+r Ts+ fT 48 


Equation (2-3) may be used to eliminate L,,(yk, s) from (2-2), yielding 


| exp {iyk-v, — k°/43} 
i ow - exp fik-v,} + oP ee — =. 48} 
s+ y T (s+r ) 


(2-4) 


> 


So " 
l exp |-¥ h 48) L,,. (yk, 8). 


+ oe 
T 8s oa T 
Replacing k by y’k in (2-2), the resulting equation may be used to eliminate L,,(y’k, s) 
from (2-4). Continuing in this fashion yields the solution 


— (—1)" exp {iy"k-v, — (k*/48)A,} 


L.(k, 3) = > 


mer 7" (s + ra ae 
where 
1 a i2n 
A. = fi 
-¥ 


It is to be noted that the series is absolutely convergent. 
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T..(k, t) may readily be obtained from L,,(k, s) by taking the inverse Laplace trans- 
form. Using a Bromwich contour it is apparent that 


rte, ) = chef” "hae, 0) de 
—. ~ [ e*'L,,(k, s) ds 


; _ 
= Residus ats = - zs ol L (k, S) 


Pe »B i exp ;24 ot = (k 138)A } (2-5) 
a S 
| tesidue at s = —-— of —_ “se | 
' erry 
Pe. “ie, = exp {iy'k vo — (k/4B)A,} ed | - ft 
ee Fo, 2 at exp ) 
7] 
1 /t t | 
-|> 4 exp |7y'k-v i 19)d,) | exp {— 
7! ¥y | T) 
The inverse Fourier transform may now be performed and this yields 
j es 
wiv, t) = B73 | exp {—7k-v}7',(k, ¢) dk 
exp {—t/7 I ) | I 
: : | exp §2(y"Vv) — v)-k - A, ? dk 
2r 7 7 1G 
J +) : | ( , / 8 \ 2 B ; a ; 
= exp = v— vo) + p Died bal | ber ] exp )— A (V — ¥°Vo) | 2-6) 
ies 1 . 7 TO \ _ ) 


This solution may be verified by substitution. 
Equation (11) may also be solved in the following way [4]. Consider the sequence 


of equations, 


ou V — WW, V / 
Ol Vv 1 v. £ ’ . 
+ | w,(v’, HAV vv’) dv 
Ey ; a 
2-7 
} V witv, / , F 
; | Vv L)Q(V yV ) dv ete 
subject to the initial conditions 
yu v, 0 d(v — Vv 
) 2 
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Plainly, 


wv, t) = >> w,(v, 2) (2-9) 


0 


satisfies Eq. (11) and the condition w(v, 0) = 6(v — vp). 
Then 
wo(v, ~) = d(v — vo) exp {—2Z/7} 


and 


/) at . ~ 
w,(V, t) exp Bz f | exp { \ | w,-1(v’’, s)\Q(v — yv’’) dv’”’ ds (2-10) 
T 0 T © 


satisfy the equations and one need only evaluate the sequence of functions, w,(v, f). It 
is seen from this last equation that if w,_,(v, é) is a product of a function of v and a 
function of ¢, w,(v, £) is also such a product. Since w, has such a form, all our w, de- 
compose in this way. Assume w,(v, ?) = U,,(v)g,(t). Then 


t) PF s 
qg,(t) = exp t.. i | exp {hon ds 


and 
U,(v) = Qo | U’,-.(v’) exp {-—B(v — yv’)"} dv’. (2-11) 
It is now easy to see that 
: ) 
t f t\ 
g(t) = — exp) ——/. (2-12) 
nN: T) 
U,(v) has the form a, exp {—8,(v — 5,)°}, and a, , 8, , 5, are connected by recursion 
relations derived from Eq. (2-11), stating 
Q exp {—£ — §.,)°} = @ | a, exp {—8,(v’ — 6,)°} exp {-—B(v — yv’)*} dv’. 
Tr} f1Ves 
Xn+1 ——; Qa, with a, =@G 
Pn + YB 
., = 3B 
J) “ (2-13) 
é = 
o — yo 
ri) "Vo 
y"—-1 
where A, = ; 
: ot (2-14) 


a= (4) 


When these are substituted into the series of Eq. (2-9), the solution is again obtained. 
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Introduction. Although several problems of heat flow in composite cylinders have 
been studied, all the cases considered treat the heat flow in the radial direction only 
[1, 2, 3]. The case of combined radial and axial heat flow in composite cylinders presents 
an interesting boundary value problem which has also considerable significance in the 
theory of vibrations and propagation of electromagnetic waves [4, 5, 6]. In this paper, 
we consider a case of combined radial and axial heat flow in the unsteady state in finite 
cylinders composed of two coaxial parts of different materials. The temperature distri- 
bution in the cylinder at any instant under the assumed boundary and initial conditions 
has been obtained by the use of the Laplace transformation. The procedure is illustrated 
by a numerical calculation in a particular case. 

The Problem. Composite cylinder made of two different materials, the inner cylinder 
0 <r < aand the outer cylinder a < r < b having thermal conductivity and diffusivity 
coefficients e, and k, and e« and k, respectively.t Boundary conditions: The flat ends 
of the cylinder x = 0 and x = l kept at zero temperature with the outer surface insulated 
and perfect thermal contact at r = a between the two coaxial parts. Initially the cylinder 
is assumed to be heated to constant unit temperature. Required the temperature dis- 
tribution in the cylinder for any time ¢ > 0 (see Fig. 1). 
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Fia. 1. Composite cylinder, the inner cylinder having conductivity and diffusivity coefficients «, and k! 
respectively and the outer e and ky. 
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Method of solution. The boundary value problem for the temperatures u,; and uz 
in the inner and outer cylinders may be stated as follows: 


ou Oru l Ou Ou 
us = 4,(2 7 = + 4) O0O<r<a, (1) 
ot or r or Ox 
Ou Our 1 Ou, OU» = 
—2 = p(s +-—+ 3), as 7T = &, (2) 
ot or r or OX / 
with the initial conditions 
uy; &, O) = wel, 2, = I, (3) 
and the boundary conditions 
u(r, 0, ) = u(r, 0, 0 = 0 (4) 
u(r, l, ) = ur, l, ) =0 (5) 
uta. ¢. 0) = 24,2, 0 (6) 
€, 0u,(a, x, t)/Or = € du.(a, x, t)/dr, (7) 


é, 0u.(6; 2, t)/0r = 0. 


We require the temperature functions u,(r, x, t) and u,(r, 2, é) satisfying the equations 
from (1) to (8) inclusive. 
Let 


U(r, z, 8) = | u(r, x, te ** dt 
be the Laplace transforms of u,; and u, . The transforms U,(r, 2, s) and U,2(r, z, s) will 


then satisfy the equations 


oU,;: 1 OU, o°U,; _ l 
aU, , 1s a Se) (9 
or” + r or + Ox ' Tae k,’ } 
ru... 0, , Us #8 « | 
——S 24 : U, = -~, (10) 
or” + rar wee Re : 
and the boundary conditions 
U,(r, 0, s) = U(r, 0, s) = 0, (11) 
Ustr; f = Ur. 1. 3) =0, (12) 
U,(a, x, 8) = Ufa, z, 8), (13) 
e, 0U,(a, x, s)/dr = €, 0U,(a, x, s)/dr, (14) 
0U,(b, x, s)/ar = 0. (15) 
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In order that U,(r, x, s) and U,(r, x, s) may vanish at x = 0 and xz = 1 as required 
by the boundary conditions (11) and (12), we expand U, and U; as well as the constants 
—1/k, and —1/k, on the right hand side of equations (9) and (10) in Fourier sine series 


U(r, z, 8) = > V,,(7, 8) sin (nrx/l), n= 13,5, --*, 


U.(r, x, 8) = } V.,.(7, 8) sin (nra/l), n= 1,3, 5,---, 
l ey 
=“; = ¥ b,, sin (nrx/l), n = 1, 3, 5, 
“1 n 


l eee 
~~ = x b., sin (nrx/l), n= 1,3,5,+->. 
For the sake of brevity, let 
a = s/k, + n'a’ /P (16) 
32 = s/k. + nx’ /P (17) 


The radial functions V,, and V»2, now satisfy respectively the equations 


PV yn 4 1 OVie _ cay, + 9s) = 0, 


or r or 


av, 1 0V, ( bss 
ore 4 = _ gt y,, + 3) = 0. 
a + a BN Von + 3 
These have the solutions 
‘ Din 
V,,(r, 8) = A,lo(a,r) — ats, 
ay 


I] 


B,Io(B.r) + C.Ko(Brr) — =, 


where J, and K, are respectively the modified Bessel functions of the first and the second 
kinds of zeroth order. The constants A, , B, , etc. are now to be determined from the 
boundary conditions (13), (14) and (15). Thus we obtain 


bin bon \ {1,(8,a)K,(8,6) — I,(8,b)K,(8,a) } Zo(an7) 
> | eB. 3 Cee 


V.,(r, 8) 








l (7, 2,83) = a in; 8. 7 A,,(s) 
(18) 
bs . NTx 
eine Fh ae 
a, l 
Ui aaw F (Bs o bas {To(Bnt)Ki(Bnb) — Ii(Bub)Ko(B,r) | Ti(.a) 
r,z#4s = €;A, . 8: A,(s) 
(19) 


bs sie nre 
- —, 
Bn l 
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oo 


where 
A, (s) = €:8,Jo(a,a)[I,(8,a)K,(8,b) — I,(8,b)K,(6,a)] 
(20) 
— ¢,a,1,(a,a)[Io(8,a)K,(6,b) + 1,(8,b)Ko(6,a) |. 


The temperature distribution functions u,(r, x, t) and u(r, xz, t) may be now obtained 
from U,(r, z, s) and U,(r, x, s) by the inversion integrals [7] 


Pier a 
u(r, 2, -) = 5— | U(r, z, se’ ds (21) 
Oat Scie 
FS oer’? ry ia 
u(r, z, t) = = U(r, x, sje" ds (22) 
2ri J - 


The integrals in (21) and (22) may be expressed as 277 times the sum of the residues 
of the corresponding integrands at their poles. In evaluating the residue of U(r, 2, s) 
exp (st) it will be noted that the first term in U,(r, 2, s) has got the factor a,8,A,(s) in 
the denominator. It will be seen further that 6, = 0 does not give rise to a pole since 
the expression remains finite (on account of the singularity of K, at the origin) as 
8, — 0 so that the only poles are those due to a, = 0 and A, = 0. However, the residue 
of the first term at a, = 0 cancels with that due to the second term and hence the only 
significant poles are those of A, = 0. Similar remarks apply in evaluating the residue 


of U(r, x, s) exp (st). One obtains therefore 


; dren ki; — kp. nme 
u(r, x, t) = y 2 cr nsin i 


(23) 
= {7,(8,,a)K,(8,;b) — 1,(B,;b)K,(B,;@)} Dolan) x,5¢ 
nn oat | “erence Sie 
On iBnj An(Ani) 
dre, k, — ky - Nx 
u(r, x, t) = — ———— ae 
7 »» r 6k l 
(24) 
5 Lola) Ks (Bai) — (usb) KulBait) TiC) 20 
, Qn jBnj An(Ani) 
where X,,; are the zeroes of A,(s) = 0, (n = 1, 3, 5, ---) and a,; and £,; are the corre- 
sponding values defined by the relations (16) and (17) when s has the values \,,; . 
The zeroes of A, may be obtained from equation (20) by solving the equation 
6a, Ii(asa) _ 11(8,0)K (8,0) — 11(8.b)K (8,0) iit 


€8, Io(B.a)  1o(8,a)K,(B,.b) + 1,(8,b)Ko(Bn,a) 


graphically. The equation (25) may be transformed into a form more suitable for 


numerical work by the substitutions 
a, = 12, B.a = ty, 8,b = ipy, (p = b/a), 
x and y being related on account of (16) and (17) by the equation 


yy = a (x? + n'a’ /l’), (26) 
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* = k,/k, and & = (k, — k,.)/k, are dimensionless constants. Introducing an- 


other dimensionless constant o”” = ¢,/¢€. and transforming the J and K functions into 


the corresponding J and Y functions* by the relations 


I,(iz) = 7 J,(2), 


where o 


K,(iz) = 7" [—J,@) + tY,(2)]x/2, 


we obtain in place of (25) 
J (x) J (py) Vily) — Jily)¥ Y, (oy) | (27) 


* Joa)” Ti(py) Vol 





Ji(py) Yo(y) — Joly) ¥s(py) 


Equation (27) has real roots and may be solved by plotting the right and left hand sides 
as functions of x. (Note that y on the right hand side is not the corresponding ordinate, 


but determined by (26)). 


Let 
e= €., j3=1,2,3,---; n = 1,3, 5, -+- 
be the roots of equation (27), the double subscript indicating that &; is the jth root 
of A, = 0. Let the corresponding values of a,,a, etc. be 
a7 = t.;, Bai@ = tom; ; Briz0 = tponn; ; 
where 


= ¢°, + (nrai/l)’ by equation (26). Then 


, = —(8,/a? + n’x’/P)k, 


u,(r, x, t) and u(r, x, t) can be now expressed in terms of &,; and 7,; as follows: 





8 2 k sae mi 9 ni Jo nj 
u(r, z, 2) = “2 ‘ ~ ~s nsin —~ > omer rw? ad r/a) gent (28) 





Sra® €; ik 3 ie ae * - sin ™ —> i 10( PT Nn | » TNniT a) J ,(&,;) errit (29) 


UNT, X, i) = /? » F £5; aj Dn(Ani) 
where 
€2 20 Nn s2 eu 
Dix. = (2: — a Gud Fy, (pom; ’ TNni) 
o oni nj 
il ni ’ 
+ ae & J (Eni) Fo0( pon: ’ on ;) 
F ni 
(30) 
+ €20d o(Enj)F 10(oMn; ’ PoNn;) 
+ i (1 saat )Juléa)Fos(one ’ Poni); 
o oC 
Fy, = J,(2)¥(y) — JW)Y,@). (31) 


*Ref. [3], Appendix III. 
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Verification of the solution.* As the series solution established by the Laplace trans- 
form method is purely formal, it is necessary to show that it satisfies all the conditions 
of the boundary value problem and is unique. It is obvious that the series solutions (28) 
and (29) for u,(r, 2, 4) and u.(r, 2, é) respectively satisfy the boundary conditions (4) 
) and (8) are satisfied by direct 


and (5). It is also seen that the boundary conditions (7 
substitution of the expressions for u, and wu, , and (6) is satisfied on account of the 
relation A,(A,;) = 0. It only remains, therefore, to verify that the initial condition, viz., 
u,; = u. = 1 for ¢ = 0. This is done most conveniently with the contour integral form 
of the solutions (21) and (22). Consider u, for example. 

For t = 0 we have from (18) and (21) 


4 . nrx ? mer 2 
u(r, z, 0) = =o = a oe. me ae ee A,1,(a,r) ds (32) 
py» nT l »» nr b Set J onic 
where 
e(k, — k.) (nr\’ 1,(B,a)K,(6,b) — 1,(8,b)K,(6,a) pa 
A, — ~ 5 2 (33) 
Relts l 8,0¢,A,,(8) 
since 
1 . naz ; 
7: sin ——- = 1 we may write u=1+y, 
7 l 
where 
4 , ee i jr" ', 
n= sin] A, Io(a,r) ds. (34) 
~ nT ( Betd oxic 
It thus suffices to show that v, = 0. The path of integration is a straight line parallel to 


the imaginary axis such that all the poles of the integrand lie to the left of this line. 
As the poles X,,; are all negative we can choose the path with any y > 0. We shall choose 
y large and positive. Since a2 = s/k, + n’x’/I’ and B% = s/k. + n’x’/P’ it is clear that 
| a, | and | 8; | will be large both for large £ and large 7. Further, if k, > k2 (say) we have 
| n . Mn Als alee oc als 7). ’ V1 v2 (Sa! ‘ 
on the path of integration | a, | < | B, 

Replacing now the modified Bessel functions in equation (33) for A, by their asymp- 
totic expansions for large argument and retaining only the dominant terms in the numera- 


tor and denominator we find that the integral in (34) becomes, apart from a constant 


(2) e ihe ” aia ds 
n\— | Pine eg 372 
r Jy-io €B, + €,a, (a,8,) 


It may be shown that the absolute value of this expression is less than 


factor 


‘i (2) “a J ( Y 4 ae] | [ c dn 

“\€; T €o) > eo) = = 4 a, Ra = 7, D2 3 3\1/3 | 372 2 

Te P\\aK, 7 21 Tl, P+ Oe)” lal 1B.) 
Thus the absolute values of the terms of the series in (34) are majorized by 

fs) : J Fr ( \ +» K( ) f nr” ( ) (35 

¢,\- exp 4 -=— (a — r) K(n) exp 4-3 (a — 35 

7) &P \—ar, Cee ae 


*Ref. [3], Appendix I. 
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where K(n) is O(nu) with a fixed finite 4. The expression (35) shows that v, can be made 
arbitrarily small by making y large. Hence it follows that v, = 0. Similar reasoning shows 
that us also satisfies the initial condition. 

The proof of the uniqueness of the solution is well known and need not be repeated 


here. 
A numerical example. To illustrate the numerical procedure, the temperature dis- 
tribution in a composite cylinder having the following parameters is calculated 


k, = j. k, = 0.1, e, = 0.4, e. = 0.04 


a = I, b = 1.5, | = 90; u(r, xz, 0) = u(r, z, 0) = 1 


Equation (27) leading to the roots &,; in this particular case is plotted in Fig. 2. 
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Fic. 2. Graphical solution of the first few roots of the eq. 27 for n = 1, 3, 5. 


The first few roots are given in the table below. 


! 2 3 | 1 
n | 
! 1.18 2.83 3.88 5.04 
3 1.03 2.74 3.85 1.98 
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Fic. 3. Distribution of temperature in the composite cylinder along the radius at x = 1/4 and x = 21/3 


fort = l andt = 2 secs. 
The distribution of temperature along the radius at x = 21/3 and 1/4 for ¢ = 1 sec. and 
2 secs. is plotted in Fig. 3. 
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—NOTES— 
ON THE RATE OF CONVERGENCE OF RELAXATION METHODS* 


By R. Piunxerr (The Rice Institute)** 


A recent paper by Frankel [1] gives the rates of convergence and a time estimate for 
the solution of finite approximations to Poisson’s equations and the biharmonic equation 
by some of the standard iteration methods. This shows that in general the time required 
is prohibitive for a reasonably large number of points. It is well known that relaxation 
methods [2] are faster for hand computing but an estimate of rate of convergence is of 
interest for machine programming for the prospective use of digital calculators. In 
general the finite approximation to partial differential equations may be written: 


Az+f=0 (1) 


where z is an n-element column matrix of the unknown values, and n is the number of 
points taken in the region. A is an n X n square matrix, with a high degree of regularity; 
the elements of the main diagonal are the largest and most of the other elements are zero; 
f isa known column matrix. If z,, is the mth approximation to z, the matrix of the residuals 
is defined by 


Az, + f = RK, (2) 


Then one element of z,, is adjusted in such a way as to make the greatest reduction in 
the value of R,, . Thus, if e; is a unit vector of A-space, i.e., e; has zero for all its elements 
except the 7-th element which is 1, 


Zm+1 = Sm — Cm4i€i (3) 


So far the value of F2,, has not been defined, but it is clear that if all the elements of R,, 
are reduced to zero z,, = z. The usual approach in the numerical solution by relaxation 
methods is intuitive, but it can easily be seen that if 


CR 
Cas = Ae, 
then 

ent. = G; 


Thus, this has the effect of reducing one element of F,,,, to zero even though it may 
actually increase the others. It has been shown [3], if A = A’ and is a positive definite 
form, that 


H,, = 32,Az, + 2f (5) 


will be reduced by such a step and the process must ultimately converge. It is difficult to 
get an estimate for H,, , however, and so this is not the easiest criterion to use. In our case 


we shall use the more customary standard of 


|R. | = (RER.)'”. (6) 


*Received February 22, 1951. 
**Now with General Engineering Laboratory, General Electric Company, Schenectady, N. Y. 
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From Eqs. (3) and (2) 
Rass = Re — Casi Ae; (7) 
and 
Pinca ee, Qem+ Ae; + Cn.1e(A’ Ae; (8) 
with no restrictions on A. Thus 
R,, | — | Runs |)(| Re | + | Ror |) = Weill Ae; — 641e5A’ Ae; . (9) 
Since in general | R,, | and | R,,,, | differ very little, letting | R,. | — | Ras. | = ARns: 


Qem+ hh, Ae, 


> 2 zee 
ARn+ = 2|R 


churetA’ Ae, 


(10) 


Let the elements of R,, be random variables with a mean value of zero, a maximum 
value of p,, and a standard deviation of bp,, where b is a number less than one. Then the 
mathematical expectation of | R,, | is: 

R,, | = n’’*bp,, . (11) 
The value of c,,.; is found from Eq. (4) by taking 7 such as to make c,,,, take its largest 
value. From Eq. (10) the mathematical expectation of AR,,,,/| R, | can be found. 
From Eq. (4), 
Cm _ Pm 


Cy = 4 = 


+1 
e; Ae, a; 


and 


Ri Ae; = Ajipm - 


where a;, is a typical element of the main diagonal. 


Likewise 


eA’ Ae, = la;; , 
where / is a small number greater than one [for Poisson’s equation 1 = 1.25, for the bi- 
harmonic / = 1.72]. 
Thus, 


Wms Ae; = 2p2 , 


2 


2 2 
Cn+ie; A’ Ae, lon 


ll 


or, the mathematical expectation of 


AR. pmiz — l 1—1/2 I 
net ae =k-, (12) 


| 22d" om 7 bn tet 
where k is a number not much greater than one unless 6, the standard deviation ratio, is 
very small. Thus for any reasonable distribution of the elements of R,, , the probability 
of AR,,.,/| R,. | being very much greater than k/n is small, or the standard deviation of 


Um 
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this quantity must be of the same order of magnitude as the quantity itself. Since because 
of the way it is found the correlation among these ratios for different m is small, then 


Ree | i( 7 Altes: | 
i TT Ra... | (13) 
and 


Rese Stag (y — Mein) ye B_1(E 
log R = » log | 1 R... a= > 9 \n 


ty» 


Thus the mathematical expectation of 


log Remsn = —k+ o(4) 
R.,, n 


and a good working value for 
> 
r= a ee (14) 


for large n. 
The ratio of the standard deviation of r, to r, is of the order of magnitude of 1/n* 
times the ratio of the standard deviation of 1 — AR,,,,/| R,, 
Since 
icon ARn+: _ i 
ei , 
and the standard deviation has already been seen to be of the order of k/n, this ratio will 
also be of the order of k/n. Thus the standard deviation of r, , which is even smaller by a 
factor of 1/n*, would introduce little chance that r, can be greatly in error. 

Since k is a number close to one, something like 20 such steps as are indicated in 
equation (14) will reduce the error as measured by | R,, | to at most 10~° of its original 
value. This is in general agreement with the usual experience in relaxation techniques [2]. 
Each such series of n steps requires about as much computation as one iteration of the 
whole matrix if the process of selection of c,,,, is ignored; this is of the order of dn arith- 
metical manipulations where d is a small whole number of the order of 10, depending on 
the number of terms of one row of A which differ from zero. If the operation of ‘‘ccompare”’ 
which is necessary to find a maximum value for c,, , takes e arithmetical operations, a 
total of about 20 (e + d)n operations will be necessary to reduce | R,, | to 10~° of its 
original value. The presently contemplated methods for “compare”’ [4] would make e 
equal to n which makes the total about 20 n° for large n regardless of the nature of A. 
This may be compared with the best results obtained by Frankel for iteration methods of 
about 20 n*”” for Poisson’s equation and about 20 n’ for the biharmonic equation. Thus, 
there is no saving in a relaxation method, which is actually more difficult to program, 
unless the problem is more complicated than the biharmonic equation. However, there 
is a possibility of making the “compare”’ operation more rapid by proper design, which 
would reduce the above estimate. From a purely heuristic viewpoint it is clear that by 
any method, at least h m operations are necessary to evaluate n points, where h is a small 
number greater than one; the difference between this and the above estimate is com- 
pletely taken up by the “‘compare’’ operation. 
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THE STABILITY EQUATION WITH PERIODIC COEFFICIENTS* 
By HIRSH COHEN (Haifa, Israel) 


In a large number of physical problems involving periodic motion, dynamic stability 
considerations result in stability differential equations which have periodic coefficients. 
In particular, if the physical system is described by a non-linear second order ordinary 
differential equation, a second order equation of the Floquet type appears. That this 
is not an isolated case becomes apparent if one reviews the large volume of non-linear 
mechanics literature of the past few years. The problem to be discussed in this note is 
even more specialized than the one just introduced but the same review through the 
literature will reveal that it is an important case. This is the stability problem which 
results when the non-linear element has small effect on the system, i.e., when the re- 
sultant motion is near to the motion of the linearized system. 

As an example consider the van der Pol equation 

y’—cdl—y)y’ +y=0 (1) 
where primes refer to differentiation with respect to ¢. If y is taken to be of 0 (1) then 
e < 1. The usual stability considerations involve the addition of a small (of order e) 
time-dependent function, »(¢), to an exact or approximate solution y)(¢). On substitution 
into (1) of y = y + v(t), the equation of first order in v(é) is 


v’’ — 1 — yo’ + (1 + eyoyiv = 0. (2) 


If the solution is to be a periodic approximation to y, then yp is periodic and (2) repre- 
sents an example of the general equation dealt with herein, namely 


ws) 


[ l 
ul? + ep(t u’ + él q(t) + u = U, (, 
€/ 


where u is the disturbance function being used to ‘test’? some system, and p(/) and 
q(t) are periodic functions of period 27/w. 

It can be seen immediately that the Mathieu equation is a special case of (3). Further- 
more, it would appear useful to remove the first order term in (3) and thus reduce it 


to at worst a Hill equation. This may be done by the substitution 


u = v(t) exp = p(t) dt }. (4 


“ 


*Received December 19, 1951. 
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with the resulting equation in v 
v’ +v+ eg,(dv = 0, (5) 


and is in fact the method employed by McLachlan [1].* 

The idea here, however, is to work directly with (3). It was found in dealing with 
the stability of subharmoniecs of the forced van der Pol equation [2] that the transforma- 
tion (4) and the resulting Hill equation (5) were cumbersome to work with when the 
desired information was only whether u(/) was a function which increased, decreased, 
or remained periodic with time. (This, of course, is what is meant here by stability. If 
the small disturbance u(/) grows with time, the original physical system is said to be 
unstable 

The following approach will be adopted: An analysis due originally to Poincare but 
used in the form given by Friedrichs [3] will be applied to (3) to discover if periodic 
analytic solutions to (3) exist. From the general Floquet theory [4] and the theorems 
of Haupt [5] for the Hill equation we are led to expect that these periodic solutions will 
form the boundaries between the stable and unstable regions. Once assured that there 
are periodic analytic solutions (analyticity in ¢ is implied by the general existence 
theorems; it is only the periodicity of these analytic solutions that is tested), the solution 
u(t) is expanded in a power series in ¢. This last named step will again produce a purely 
periodic solution but will also produce the conditions on p(t), gi(Q and for which stable 
solutions exist. The feature of this analysis which may be novel is that it is shown that 
the periodic coefficients need not have exactly the linearized period of (3) but may be 
somewhat different [according to the form of p(t) and q,(Q] and still produce a stable 
solution. 

It should be emphasized here that this is intended to provide a quantitative study 
of the special equation (3) dealt with, and even that only in the restricted region of « 
small. Qualitative investigations begin with Liapounoff [6] and have been taken up by 
other authors. 

In order to use Friedrich’s approach directly, let ¢ = wf and consider a phase shift 
6 in o such that the equation (3) now written in the form 


u’ tu = —elp(O+ du’ + g,(6 + Su] 
has the special initial conditions 
u(0) = A, u’(0) = 0. 


Following Friedrichs we introduce the variable such that the period of the resultant 
motion, 7, equals 2x + en(e). Friedrichs then seeks values of A, 6, and 7, taken to be 
analytic in e, such that periodicity conditions on the solutions near « = 0 are satisfied. 
The method is given fully in the textbook of Stoker [8, p. 233] and need not be repeated 
here. In the ease of a linear equation, A is independent of wf and is not involved in the 
periodicity considerations. The condition for periodicity of solutions finally obtained is: 


\ 
1/2) 


No = oe a k (a. — d,)” 4 (Co + b.)° — a * 


*Numbers in brackets refer to references given at the end of the paper. 
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a2 ] 


Here 7 is the first term in an expansion of n(e) in powers of ¢, and the ¢ , % , 
b. , C2 , d> are the coefficients in the Fourier series expansion of p and q, , namely 


p(t) = a + 7 (a, sin nwt + b, cos nwt), 
i 


q(t) = Co + >> (c, sin nwt + d,, cos nwt). 


— 
l 


Translating back into terms of w, we have 
en . f l l ont 

i / — 2 2 2 \ 

w=l1- =1—-— jo += (a, — d.) + (co + b.)° — ag ?. 


).. ¢ 
- 2 | J 


Thus periodic solutions will exist when the period of the coefficients is given by 


i? 


a / ( 1/2) 
: T=-—= Onl 1 a < $C) + E (a2 — do)” + : (Co + b,)” = a | ) 
w \ \ i f J 


Let us now investigate (3) by the customary formal expansion in a power series in e. 


With the change of variable r = wt, we have 
? € , | € ! I . 
NT) Pr DTT) or Aaa) u(r) = O. (6) 
WwW € 


It is well known [4] that the fundamental solution to the Floquet equation is given by 
u = e'g(7), 

where X is a complex number and g(r) has period 27. Let us take A, g, and w to be analytic 
functions of the small parameter e so that we may write the expansions 

A=rA tA téerXx,+°°:,7z 

9=PJtentege + 

wo= 1 +e, + ew + 

Substituting these into (6), equations of any desired degree in ¢ are obtained. For 

the degree zero: 

gs’ + 2rogs + Ac + I) = O. 
Now if the condition that g be periodic is imposed, the g, must be periodic; in particular, 
go must have the linearized period, 27. This will obtain if 

MY=+tt, g=cte””. 


Here one of the integration constants has been taken to be unity with no loss in gene- 
rality. Proceeding, the equation of first degree in e can now be found making use of the 
zeroth order solutions, A) and go . The result is 

; -2roT ¢ ¢ -Ao2r — 

gi’ + 2rogi = —(c + ¢ *°T\(QiAo — Zw. + Gi + Aop] + Woe (2A, + pl]. (7) 
Since the solution g, is to be periodic, the terms on the right hand side of (7) which 
will produce solutions that are non-periodic must be eliminated. This represents the 
imposition of the periodicity condition as a boundary condition. The method, often 
termed the casting out of secular terms, was given by Poincare and used extensively by 
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Duffing and many others. The terms on the right hand side of (7) which give rise to 
secular solutions are the constant terms and the terms with e~”*’ as multiplier. Again 


we use an expansion of the periodic functions: 


| 
4 
| 


=a, + >> a,sin nr + b, cos nz, 
n=1 


g(r) = Co + > De c, sin nr + d,, cos nr. 
n=1 


‘ . » 2 
Collecting constants and coefficients of e“**’, we have 
Co AoA d, Aobe 
— =" i _— 2. & 3 
c[2A1Ao Qu, + Cy + ApAo] ry + 95 9 +. . = 0, 





Co d, Ande dodo 0. 


eT an — — . at mal aie an ’ = 
[ 2ZdoAj AoA aw) + Co Al 21 + 2 23 + D 


Eliminating c we obtain 


2, = —a+ k (ec. + bo)? + ; (a, — d.)? — (cy) — 2w;)’ : (8) 





To check the previous result, we set 4, = O and obtain 


. 1 s 1 . P 1/2 
W; ~ +; E (Co + bo)” + 4 (a, — ad.) — a | 3 


~ 


Since w = 1 + ew, in this approximation and w = 1 —em)/2zm has been used, then w, = 


—n/2zr. The result is 


— —r(c, + E fe + hy + } (a, — d,)* — a ). 


Returning to (8), it is observed: 
1) if 
2 l 2 2 
(Co — 2w,)° > { [(C2 + be)” + (a2 — d.) |, 
stability is determined entirely by ao 
a) do < 0, A, > O, and the resulting w(¢) increases with time 


and u(t) decreases with time. 


db a>ba, < 9, 


(cy — 2a)? < F {le + bs)* + (@ — d,)*}, 


stability is determined by 


> 


— os 4. E {(c, + b.)? + (a, — d.)*} — ( —- 2a) ; 


“< 


Let us consider as an example to explicate the above work the equation (2) where 
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Yo Will be taken to be A sin ¢. (This is exactly the problem solved by McLachlan in [1], 


p. 190). Then 
( 3 7 ‘7 
[ 9 — I — 9 cos 2, 


Pb) = 
qth = A’® sin 2t, 
A* 
&®= =>! 
a = ¢ = d, = 0, c, = A? b, = —4 


Notice that for periodicity we would need 


w =i(4 _A s#-a) <6 


2 \16 ! 


since here w = 1. If A*® = 4, then this condition is satisfied. Further, using (8) since 
Co — 2w, = Oand 3{ (co + b.)” + (a, — d.)"} > 0 
A* ] 
mee Ts 


Thus \, < 0 for A® > 4 and >0 for A~ < 4. 
This example as was remarked is given by McLachlan [1] and was presented merely 
to show the ease in which stability characteristics may be obtained once A, is computed 


in terms of dy), 2, bo , Co, and d2. 


REFERENCES 
1. N. McLachlan, Ordinary non-linear differential equations, Oxford, 1950. 
2. H. G. Cohen, Subharmonic synchronization of the forced Van der Pol equation (To appear in the Pro- 
ceedings of the Colloquium on non-linear vibrations, Ile de Porquerolles, August, 1951. 
3. J. J. Stoker, Non-linear vibrations, Interscience, New York, 1950. (see especially Appendix I). 
4, E. L. Ince, Ordinary differential equations, Dover, New York. 
5. M. J. O. Strutt, Lamésche, Mathieusche und verwandte Funktionen tn Phystk und Technik, Julius 
Springer, Berlin, 1932 
A. Liapounoff, Probléme general de la stabilité du mouvement, Princeton University Press, Princeton, 


1949. 


ON THE RELATIONSHIP BETWEEN THE MARTIENSSON AND DUFFING 
METHODS FOR NONLINEAR VIBRATIONS* 


3y ROBERT E. ROBERSON (Mechanics Division, Naval Research Laboratory) 


The background for a number of one-term approximation methods and their appli- 

- . . . * ‘ ° 1 
cation to forced nonlinear vibrations has recently been discussed by Schwesinger. 

*Received Aug. 15, 1951. This paper corresponds to part of a dissertation submitted to Washington 
University in partial fulfillment of the requirements for the degree of Doctor of Philosophy. 

1G, Schwesinger, On one-term approximations of forced nonharmonic vibrations, J. Appl. Mech. 17, 


202-208 (1950). Note that he attributes to Riidenberg the method that is designated here as Martiensson’s 


method. 
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As he points out, such one-term approximations may be useful in problems of analysis 
when the response is close to sinusoidal, as is frequently the case for small nonlinearities. 
They are not the ultimate tool, however, for one can always go to approximations which 
contain the higher harmonic terms as well. The situation is different for problems of 
nonlinear synthesis. There, the nature of the problem is such that approximations of 
more than one term are virtually prohibited, and one is forced to accept answers that 
are given by the one-term methods. Thus, there is a body of problems for which such 
methods have intrinsic importance. 

The requirements of synthesis have led to a re-examination of some of these one- 
term methods. For this purpose, it is important to know in a general way whether the 
use of the various methods leads to different synthesis. In particular, one should know 
which of the methods are really independent of one another, in order to know how many 
synthesis possibilities must be examined. It was found that the results obtained by the 
methods known by the names of Martiensson (or Riidenberg) and Duffing’ are not 
independent. For the class of equations to which both methods are applicable, there is 
a simple relationship between their results. The purpose of this note is to develop this 
relati mship. 

The fact that the results of the two methods are simply related is the more striking 
because their rationales are so different. Briefly, these are as follows. In the Martiensson 
method, since an assumed one-term solution does not generally satisfy the differential 
equation identically, one forces satisfaction at two points, say ¢ = O and at the quarter- 
period. In this way, an algebraic equation is obtained for the amplitude of the one- 
term solution. To use the Duffing method, one writes the equation in a form for iteration, 
and puts an assumed one-term solution into this equation as a first approximation. In 
order to enforce the requirement that the next approximation be periodic, one is obliged 
to choose a certain coefficient equal to zero. This again leads to an algebraic equation 
for the amplitude. As a development of the Lindstedt perturbation method, the Duffing 
method is perhaps the more rational procedure of the two. (It is fundamentally different, 
of course, in that it permits extension to higher approximations.) 

Now let us consider the methods in greater detail. They are usually illustrated only 
for systems with one degree of freedom, but it is not difficult to extend them to certain 
higher order systems. The Martiensson method, as ordinarily applied, is limited to 
systems without dissipation. The most general equation to which it seems appropriate, 


at least without major changes in formalism, is 
OY 2+ Lr + vf(x) = F sin t. (1) 


Here, x(/) and f are (n X 1) column matrices, F is a constant (n X 1) matrix, and L 
is a constant (n X n) matrix. Both ’ and » are scalar parameters, the former being 
related to the frequency of excitation and the latter being the nonlinearity parameter. 
Suppose that f(x) is a continuous odd function, with f(0) = 0. We will seek a periodic 
solution with period 27. 

Let us apply the Martiensson method formally to Eq. (1). We assume an approxi- 
mation solution 2°’ = asin ¢ (a being a column matrix of amplitudes), which satisfies 

2We restrict our consideration to the first approximation by this method, in effect making it a one- 


term method. 











NOTES [Vol. X, No. 3 


272 
the equation at time ¢ 0. If we require that it also satisfy the equation at / = 7/2, 
we obtain 
(2° — Da — vf(a) + F = 0. (2) 
On the other hand, let us write Eq. (1) in a form for iteration as 
(2? +a) = (2 — Da” — vf(x"’) + Fsin t. 

Using the same approximation as before, 

(2°? + 2) = [(2? — Dat F]sin t — vf(asin 2). (3) 


Since f is a periodic function of ¢ with period 27, we can write 
f(asin t) = §,(a) sin t + > Fom-1 sin (2m — 1)t. 


Following the Duffing iteration method, we require that no secular term arise when 
Eq. (3) is solved for x‘’’. This means that we must put the coefficient of sin / in Eq. (3) 
equal to zero, namely 

(2? — Lia — (a) + F = O. (4) 
Equation (4) is the analog of Eq. (2), and is identical except that the first Fourier co- 
efficient of f(a sin /) replaces the value f(a) 

This result has limited usefulness in its general form. However, in systems which 
contain only one nonlinear element, f is particularly simple. It has the representation 
f = Mog(z;), where M is a column matrix and ¢ is a scalar function of only one of the 
z-components. As perhaps the simplest possible example, we may consider the case 
where ¢(z;) = x; , ie. a simple power function. Equations (2) and (4) become re- 


spectively 


(2? — La — va + F =0 


(2° — Lia — ca’. +F=0 
where c is the constant whose value is the first Fourier sine coefficient of sin” ¢. For 


9 


example, if m = 3, thence = 3/4. 

The result of this simple special case has a very useful implication in synthesis 
problems. It means that if a synthesis is attempted for a system containing a single 
nonlinear element which obeys a power law, the system being described by Eq. (1), 
then the procedure follows identical paths for the Martiensson and Duffing methods. 
It is known at once that the optimum nonlinearity by one of the methods, say the 
Duffing, is just 1/c times as large as those predicted by the other, and that optimum 
values of any parameters contained in LZ are identical by the two methods. Nothing 
should be inferred, of course, as to which method is better as a one-term approximation, 


since this is an entirely separate problem. 





Rae? 


? 
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NOTE ON AERODYNAMIC HEATING 
WITH A VARIABLE SURFACE TEMPERATURE* 


By A. E. BRYSON (Hughes Aircraft Company) 


Emmons [1], has considered the problem of an insulated flat plate of infinite extent 
started impulsively from rest in a viscous, incompressible fluid. One of the interesting 
results of his analysis was a simple expression for the temperature recovery factor at 
the plate surface. Another interesting result can be obtained from the same problem 
by considering, instead of an insulated plate, a plate with a surface temperature that 
is a given function of time. 

As Emmons has shown, letting y, the viscosity coefficient of the fluid, p, the fluid 
density, and k, the thermal conductivity of the fluid, be constant, it follows from the 
equations of motion and the boundary conditions that the pressure is constant and the 
velocity normal to the plate is zero. The momentum equations reduce to: 

Ou ou 

at oy’ (1) 
where u is the velocity component of the fluid parallel to the surface of the plate, ¢ is 
the time, y is the distance perpendicular to the plate surface, and vy = y/p. The energy 


equation becomes: 


PP oT a 2 ee ia] (2) 


at ‘ ay?" \dy/’ 


where T is the fluid temperature and Cp is the fluid specific heat per unit mass at con- 
stant pressure. These equations are to be solved with the following boundary and initial 


conditions: 


u(O, t) = 0, u(y, 0) = U, (3) 


T(0, t) = T.(d), T(y, 0) = T.., (4) 


where U is the free stream velocity, 7'.. the free stream temperature, and 7’, the plate 
surlace temperature. 
The problem defined by Eqs. (1) and (3) for the velocity diffusion is well-known, 


the solution being 
u = Uerfly 2(vt)' 2% (9) 


Substituting (5) into (2), we have 


oT oT U* exp (—y’/2vt) ' 
--a gee, (6) 
ot oy” rp t 
where a = k/pCp. The solution to (6) with the boundary and initial conditions (4) is 


obtained by the method of heat sources and sinks and the result is: 


*Received ( Yetober 26, 1951. 
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T —-T., = Ax” | [7.(t — y°/4aB") — T..] exp (—8°) dB 
» Fs (7) 
| wy | | +0 | 
‘ ex —— — ex _ 
° e" exp (—X°/2yz7) fa(t — 7) 4a(t — 7) 
UP fay fi esp (=n Boe) “PL dalt dat = 0) 9, 
rCp Jo Jo T [4ar(t — 7)] 


The heat flow through a unit surface of the flat plate, Q, is given by the derivative 
of (7) evaluated at the surface of the plate: 
oT 1/2777 1) 1/2 i 1/277 
Q=k > = —k(rat)’*[T.(0) — T.] — k(ra) | (¢— 7) Thr) dr 
C _ “0 
7 lene (8) 
+ k(mat)~*’*r(c) U* /2 'D, 


where the recovery factor, r(c), is given by: 


i {1 — o/2\'” F 
1/2 \tan acne? ; <2 
a/2 
) (9) 


| ( a/2 
Tio) = —— 
r\|}1—oa/2 
9\1/2 ¢ 1/2 ‘ 
log [(¢/2)’" + (¢/2 — 1)”’]; o¢>2 
and ¢ = v/a = Prandtl number. Equation (8) can be written more compactly in the 


form of a Stieltjes integral as follows: 


« 


Q = —k(rat)"'”” | (1 — 7/t)'” d[T.(7) — T.], (10) 


where 
T. = T. + r(o)U"/2Cp. (11) 
The surface-temperature variation to give a prescribed heat flow variation can be 
obtained by inverting the Abel integral Eq. (10); this gives 


T, — T.(t) = (a/n)'*k™ | (t — r) Q(z) dr. (12) 


For a constant rate of heat flow this reduces to 
T. — T(t) = 2Qk‘(at/x)'”’. (12a) 
Emmons has already given this latter solution for constant heat flow to the plate [1]. 
By replacing ¢ by z/U in (1) and (2), we have the linearized boundary-layer equations 
for steady, viscous, incompressible flow past a flat plate (sometimes called the Rayleigh 
equations). The boundary conditions (3) and (4) become 


u(0, x) = 0, u(y, 0) = U, (13) 
T(O, x) = T(z), T(y, 0) = T. (14) 
These imply a semi-infinite flat plate with an arbitrary surface temperature 7’,(z), 


where z is the distance from the leading edge in the direction of flow. Therefore, an ap- 
proximation to the heat flow to a semi-infinite flat plate in a steady flow of velocity U, 
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and with an arbitrary surface-temperature distribution 7',(x), can be obtained by re- 
placing ¢ by 2/U in Eq. (10); this gives 
Q = —nr '?k(U/v2)'"0'”? (1 — &/x)'” d{(T.(é) — T.] (15) 


« 
“0 


which, for constant surface temperature reduces to the familiar form 
Q = 9 ?kK(U/a)'?o'*(T. — T,). (15a) 


Lighthill [2] has given an expression for heat transfer to an arbitrary two-dimensional 
surface in terms of the skin friction and temperature along the surface, by using the 
Fage and Falkner linear approximation of the boundary-layer profile and neglecting the 
viscous dissipation terms in the energy equation. For the case of the flat plate, his result is: 


Q = —0.339k(U/vx)'7o"* | fl — (&/x)**)"'? d[T.() — Te]. (16) 


He argued that the effect of viscous dissipation is taken care of by replacing 7... in the 
above expression by 7’. + o}U*/2Cp, the boundary-layer equilibrium temperature. If, 
following the suggestion of Lewis and Carrier [3], we replace U by 0.35U to approximate 
a mean convective velocity in the boundary layer the constant multipliers of Eqs. (15) 
and (16) are nearly equal. 

The differential equations used here and by Emmons in [1] apply to a fictitious fluid 
of constant pressure and density, but variable temperature. The equations are really 
of interest only because the compressible fluid boundary layer equations can be reduced 
to their form by the Von Mises transformation and the assumptions that yu is proportional 
to the temperature and o is constant (see for example ref. [4]). If enthalpy is used as 
the independent variable instead of temperature, no additional assumption need be 
made on the variation of the specific heat with temperature. Then the only change in 
the previous differential equations is to replace y by n where 


oy 
7 = | £ dy 
Jo Pa 
The expressions for heat transfer rate are unchanged, although as they stand Cp must 


be assumed constant. 
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A NOTE ON THE DAMPING IN ROLL OF A CRUCIFORM WINGED BODY* 


By JOHN W. MILES** (Fulbright Lecturer, Auckland University College, Auckland, New Zealand) 
1. Introduction. In the following we apply the slender body theory of Ward’ to the 
calculation of the damping in roll of a slender body of revolution of radius a, bearing 


cruciform wings of total span 2b. 
It is required to find a solution, ¢(s, z, y), to Laplace’s equation satisfying the 


boundary conditions 


$,(8, 2, ¥) |yeo = 2, a<|z| <b, (la) 
¢.(s, 2, y) |z-0 = —DY, a<|ly|<b, (1b) 
¢,(s, 7 cos 6, rsin 6) |,-, = 0, (1¢) 


where (s, 2, y) are a set of right handed, Cartesian coordinates with s measured positive 
downstream from the body nose, and p is the angular velocity about the s axis. The 


rolling moment is then given by 
N = —pU $ o(l, x, y)(x dx + y dy), (2) 


where the integral is taken around the cross section at the trailing edge (s 1), the 


latter being assumed to be straight and transverse the free stream. Whereas in Eq. (1) 


a and 6 are functions of s, exhibiting a monotonic increase therewith, we hereafter refer 


only to their values at s = 1. 
2. Solution for potential. The conformal transformation 
(z + a’/z)? — (b + a’/b)*? = (¢ —c/)’, (1) 
c* = 1(b° + a‘/b’) (2) 


maps the profile of the cruciform winged body in the z2(=zx + zy) plane on the circle 


¢| = c in the ¢ plane, the wings appearing as four, symmetrically disposed arcs of 
subtended angle 2y, , where 


‘ 


cos 2g, = (a/c)” = 2(a/b)*[1 + (a/b)*]™. (3 


Transforming the boundary conditions (1.1) and carrying out the solution to the then 


classical problem, we obtain the potential on ¢ = ce‘® in the form 
I 
‘ 2 - "i ¢ ‘ Ss 9 1/ 
@¢=>-— 2pc 7) sin 49 | {(COS 2y — COs 2¢, )+ (cos* 2y — cos’ 20) 7 
(cos 4y — cos 4y)' dy. (4) 


Carrying out the required integrations in Eq. (4), we have 


*Received August 21, 1951. 
**On leave from University of California, Los Angeles; formerly Consultant, U. S. Naval Ordnance 


Test Station, Inyokern, California. 
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| a 

; tan (¢ 0) | | sin 2 Do 
@ = —(pe’ 25) co 2¢ In | - (e + go) | _ COS 2g In | om Be + 2 + 0) | 
| tan (p — ¢) | | sin 2(9 — @) | 


+ K(sin 2y)) sin 4g + 2[E(sin 2¢,)F(sin™ (sin 2¢/sin 2g), sin 20) 
— K(sin 2¢,)E(sin™ (sin 2¢/sin 2g»), sin 2g») ](sin’ 2g) — sin” 20)", (5) 


where K and E, with single arguments, denote complete elliptic integrals of the first 
and second kinds; respectively, and F and E, with two arguments, denote the corre- 
sponding incomplete integrals of modulus sin 2g, and amplitude sin™' (sin 2¢/sin 2¢). 
For the special case of no body (a = 0, g) = 2/4) Eq. (5) reduces to 


1 + sin 2p | (a = 0) (6) 


_ -_ b? 20S fa) a 
od pb’/2m) cos 2¢ in| 1 — sin 2¢ | 


3. Rolling moment. The integrals subsequent to the substitution of Eq. (2.5) in Eq. 
(1.2) appear to be intractable, and it is expedient to proceed by substituting instead 
Iq. (2.4), which leads to 


@ 


N/N, = (2/m)*[1 + (a/b)*? >> (2 /n), (1) 
1 

No = —xpUpb* 4, (2) 

l= | sin (4ny)d[cos 2p + (cos? 2g — cos” 2y,)'”]. (3) 


The reference moment, N, , is twice that acting on a single wing of (small) span 2), so 
that N/N, represents an overall interference factor. For values of (a/b) near unity the 
convergence of Eq. (1) is poor, but the potential may be expanded in powers of ¢ to 


obtain 
N/N, = 8{1 — (a/b)? {1 — 1.57[1 — (a/b)]} + Of{[1 — (a/d)]*} (4) 


For the special case a = 0, the substitution of Eq. (2.6) in Eq. (1.2) yields an inter- 
ference factor of (8/x°), which may be checked by summing Eq. (1), viz. 


N/No = (2/m)? DU (ae) "| aad " 


in agreement with the result obtained by Adams.” We remark that the first three terms 
in Eq. (5) yield the correct result to better than 2%, implying a satisfactorily rapid 
convergence of Eq. (1) for small (a/b). 

More complete numerical results are to be given in a forthcoming NOTS report,” 


representing work supported by the Office of Naval Research. 


REFERENCES 
(1) G. N. Ward, Supersonic flow past slender pointed bodies, Q. J. Mech. Appl. Math. 2, 75-97 (1949). 
(2) G. J. Adams, Theoretical damping in roll and rolling effectiveness of slender cruciform wings, NACA 


TN 2270 (1951). 
(3) J. W. Miles, On the damping in roll of a slender cruciform wing body, USNOTS Report, Inyokern, 


California, (1951). 








278 NOTES [Vol. X, No. 3 


THE ELASTIC AXES OF A ONE-MASS ELASTICALLY SUPPORTED SYSTEM* 
By J. J. SLADE, JR. 


When an elastically supported rigid body is subjected to the action of a rectilinear 
sinusoidal force, the resulting steady motion generally consists of rectilinear and torsional 
oscillations with frequency equal to that of the exciting force. It is desired to determine 
the location of the exciting force so that the torsional oscillations are suppressed or, at 
least, so that the amplitude of these oscillations is reduced to a minimum. The problem 
arises, for example, in connection with unbalanced machines on elastic foundations, as 
well as in investigations of the dynamic characteristics of elastically supported rigid 
assemblies by means of induced vibrations. 

The two-dimensional problem has been considered under simplifying conditions.’ 
The three-mass mechanical oscillator® that produces a force the axis of which may be 
made to coincide with any line in a fixed plane, when the oscillator is in a fixed position, 
presents problems that require an extension of existing theory. The present investigation 
deals with the general case. 

We consider a rigid body of mass m that can move freely under general linear elastic 
constraints with linear damping. Let r be the displacement of its center of gravity with 
respect to its position in static equilibrium. Since only small oscillations are considered, 
elastic and damping reactions may be taken to be fixed to a moving frame with origin 
at r. 
Let d + e&, and V + e¥, be dual dyadics such that —(@ + e&,)-r is the motor’ of 
the elastic suspension and - (¥ + e¥,)-7’ that of the damping system, due to a rectilinear 
displacement, the prime denoting differentiation with respect to time. 

Finally let f be the exciter force and p a point on its line of action. It should be noted 
that in all cases considered the exciter is rigidly connected to, and forms part of the system. 
The exciter force is strictly fixed in the moving frame. 

The motion of the center of gravity of the body is governed by the equation 


mr’ + Wer + O-r = f. (1) 


There is also the moment 
c=pXf— (r+ W%-r’) (2) 
that tends to produce torsional oscillations. 
If the angular frequency of the exciter force is w, we may write 


t 


f,r,¢ = (F, R, Ce’ 


where 


OO 


(—mwa ll + wv + %)-R=F (3) 


*Received Oct. 18, 1951. 


See, for example, E. Rausch, Machinenfundamente und andere dynamische Bauaufgaben, Ch. III, 


V.D.I., Berlin, 1936. 
2R. K. Bernhard, Study on mechanical oscillators, Proc. Am. Soc. Test. Materials 29, 1016-1036 (1949). 
3],. Brand, Vector and tensor analysis, Ch. II, J. Wiley & Sons, New York, 1947. 








1952 J.J. SLADE, JR. 279 


and 
C=p X F — (®% + toW,):(—mwa 1 + iw + 6) '-F 
(4) 
=pXF+ ("+ iA)-F. 

Our problem is to determine under what conditions, if any, we can make C = 0 
with F ¥ 0. 

An axis with which the line vector F + ep X F must coincide to satisfy these condi- 
tions fully is here called an elastic axis* of the system. An axis of fixed direction with which 
the line vector must coincide to make | C | # 0a minimum will be called a quasi-elastic 
axis. 

Suppose first that the system is conservative, so that V + e¥) = 0. If oscillations about 
the axis of the free vector a are suppressed, then a-C = 0; or, since in the conservative 
case A = Q, 

0. (5) 


II 


apXF+a-T-F 

Now, the left hand member of this equation is the moment of the fixed motor 
a + ea-T about the axis F + ep X F Whence: 

Oscillations about an axis a are suppressed when the line vector F + ep X F coincides 
with a line of the null system of the motor a + ea-T. 

Let €; , @2 , @; be unit vectors in the directions of the principal axes of the elastic 
suspension. In this presentation the diagonal elements of T are zero and e,-e,-T = 0, 
so that the motor e, + ee,-T is a line vector. We therefore have the following results. 

|. Rotational oscillations about a principal axis are suppressed when the line of 
action of the exciter force coincides with a line of the special linear line complex the axis 
of which is ey + ee,-T. 

2. Rotational oscillations about two principal axes are simultaneously suppressed 
when the exciter force coincides with a line of the linear congruence, the directrices of 
which are e; + e;-T,j = k, l. 

3. The elastic axes of the system are the lines of the regulus the directrices of which 
] 9 2 


my 


are e, + ee,°T, h 
When the system is not conservative, the following equation must be added to (5): 


a-A-F = 0. (6) 


1. Oscillations about @ are suppressed when the force coincides with a line of the 
null system of a + ea-T that is perpendicular to the fixed couple vector a- A. 

Assuming that the principal axes of V coincide with those of #, as they generally do 
in practical cases, we may further state. 

1. Rotational oscillations about a principal axis e, are suppressed when F + ep X F 
is a line of the plane through e, + ee,-T perpendicular to e,- A. 

2. When F + ep X F is the line of intersection of two such planes oscillations are 
simultaneously suppressed about the corresponding two principal axes. 

In general the non-conservative system possesses no elastic axes. Since T and A are 
constants, when w is fixed, we see from Equation 4) that, if F is held fixed, | C | is a 
minimum when p is so determined that 

pXF+Tr-F=0. (7) 


‘Rausch, loc. cit., uses the terms elastiche Hauptachse and elasticher Mittelpunkt. 
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This leads to the following result: The quasi-elastic axes of the system are the lines 
of the regulus the directrices of which are e, + ee,-T, k = 1, 2, 3. 

The traces of the directrices e, + ee,-T on the principal planes have been called 
elastic centers. The locus of an elastic center starts, with w = 0, at a point that depends 
on the parameters of the elastic suspension and ends at the center of gravity (w = ©). 

In the conservative case this locus is the outside section of a hyperbola, the inside 
section corresponding to w < 0. The locus is a 4th degree algebraic curve in the non- 
conservative case. The reduced system in which one reaction goes through the center of 
gravity and the other two lie in a plane through this center has been considered in detail.” 


HEAVY DISK SUPPORTED BY CONCENTRATED FORCES* 

By YI-YUAN YU (Washington University, St. Louis, Mo.) 
a two-dimensional light disk subjected to an 
arbitrary number of concentrated forces by means of his method of complex variable 
(1, 273-274].** When the weight of the disk has to be taken into consideration, the prob- 


“ato” 
lem may still be solved in a similar way. MuscheliSvili’s notations will be followed 
il ones will be defined as they first occur. 


MuscheliSvili solved the problem of 


throughout this paper, and only addition: 
In plane problems including body forces due to gravity, the stress function U may 

still satisfy the biharmonic equation if it is defined by the following equations: 

a°U ’ au , aU 

= — Vi, a eo Fs. a (1) 

Ox OY OX OY 


in which V, is the body force potential due to gravity and is equal to wy when gravity 
acts in the negative y-direction [3], w being the specific weight of the material of the 
body. Hence, the function U may be expressed in terms of two analytic functions as 


shown by Muscheli§vili [2, 284]. 
The boundary conditions 
dz, dy dx 


dy 


= Tr, fr ‘ 
"7? dg Ge 


7. 7 
“" ds ” ds 


however, may be shown to lead to some different result. When stress components given 
by Eqs. (1) are substituted into these conditions and computations carried out in the 
(2, 301-302], the following result is obtained: 


same manner as given by Muscheli8vili 


. 


o(2) + #92) + Vi@ =i] (re + ir) ds— | Vide 


d 


If we define 


. . 


fitif.=1 | (r. + ir,) ds — | V, dz (2) 


d e 


5R. K. Bernhard and J. J. Slade, Jr., On the elastic center of one-mass plane oscillatory systems (un- 
published). Dynamics Laboratory, Bureau of Engineering Research, Rutgers University. 
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then the boundary equation for the first fundamental boundary value problem is 
gi(z) + 29;(2) + (2) = fi + tf 


which holds true on the boundary C of the body in the original z-plane. When C is 
mapped into the unit circle y in the ¢-plane by means of the function z = w(f), the 
boundary equation becomes 

w(o) 


g(c) + g'(o) + vio) = fi + the (on y) (3) 


w (a) 


Thus, except with a different definition of f, + cf. , this equation assumes the same 
form as the one for zero body forces [2, 294]. 

Modified expressions for stress and displacement components may similarly be 
derived. Only stress components in curvilinear coordinates are given here: 


Top + Too = 2H(e) + O()] + 2V(9) } 


150 — Tp + ity» = —— [a(H)#"(s) + oooh 


i] which 


Vit) = V,@(h)) = V2) 


By comparing with MuscheliSvili’s original formulas [2, 312], it can be seen that 2V(¢) 
is the only additional term due to body forces. 

The problem here concerned is that of a heavy disk having radius R supported by 
an arbitrary number of, say n, concentrated forces (X, , Y,), (X2, Yo), --+ , (Xn, Ya) 
at points on its boundary corresponding to o, , a2, «++ , ¢, respectively on the unit circle 
as shown in the figure. Obviously the supporting forces must satisfy the following con- 
ditions: 


) ae es > Y, = wah? 


(x,,¥) A 


©) 








O;) On) 


(X2N,) 


22 
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Since the boundary of the disk is free from stress everywhere except at the supporting 


points, expression (2) becomes 


fi tif. = -| V,.az= : tw / (z — 2) dz 


For a circle of radius R, we have 


2 = a(t) = Rt 

hence, 

- . ] ° j > l ‘ > o 

fit ife = 5 Ww R*(o — a) do = = wh \ 5 — loga 
The boundary equation (3) now takes the form 

: | 
; eee i ” 
o(c) + aoo'(o) + Wo) = 5 wh a log o (on ¥) (2) 


By modifying those obtained by MuscheliSvili for the problem of a light disk [1, 27 t], 
the two analytic functions for the present problem can be written down as 


| — 
Qt) =- z. X;, + 12) log(o, —- OH +¢(¢ 


ic — | (X, + 7Y,)o; 


Dae 
“= Tl 


in which ¢°(¢) and y°(¢) are functions analytic in the entire region inside y and have 


the forms [1, 272] 


=e 
| 


The other terms in ¢(¢) and ¥(¢) account for the singularities in the solution due to con- 


centrated forces and therefore have the same forms as those for a light disk. 


Substituting expressions (6) into Eq. (5) now yields, after simplification, 
iwR” l wey, nee a 
g(o) + o¢ (oe) + vio ro. =F ‘2 X, — tY;,)o.0 (7) 


from which ¢°(¢) and W°(¢) can be determined. Following the established procedure of 
the method, we formulate an integral equation by multiplying Eq. (7) through by 
1/2mi do/(o — ¢) and integrating around y. The integrals in the equation thus obtained 
ean be evaluated by means of the theorems developed by MuscheliSvili [1, 269]. The 
result gives 


‘ : ‘ iwh* ) | * - . 
C + 2(a; - 10>) +a,— 18h =~- C } oF — 1) o.C 


4 2r 


l 


in which the constant terms may be neglected. The constant a, is determined by differ- 
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entiating the rest of the equation with respect to ¢ once and setting ¢ equal to zero. 

We have finally 

iwR™ 2 l 
; *  & 


> (X, — i¥og (8) 
k 


Substituting this back into Eq. (7) and formulating the conjugate of the result, we obtain 


és) = —+ Go” + ; > ((X, — i¥Jo. — (X + iY Da Io 
T ez} 


9 
Multiplying this through by 1/277 do/(o — £) and integrating around 7, 
y(t) = —iiwR’ (9) 


Thus the problem is completely solved. The solution consists of ¢(¢) and ¥(¢) as given 
by (6), and ¢°(¢) and y°(¢) are given respectively by (8) and (9). 

The problem of a heavy disk resting on a horizontal plane was solved by J. H. 
Michell [4] and represents a special case of our problem in which 


n= 1, X, = 0, Y, = wrk’, 1 =-i 


The analytie functions reduce to 


iwh , " . , C0 os wh? 

A) = -> log (F + 1) FG > s 

wh l twh* , j ° - »2 

oe ees log ( + t) — 5 wh 
The sum of the normal stress components can be computed according to the first of 
Eqs. (4); thus, 

Top T Tee = wh —-2- = r z ) 
\ p + 2psin 6+ 1 


It can readily be shown that both normal stress components vanish at all points on the 


boundary of the disk except the point of support. 
The problem of a heavy disk resting on the ends of its horizontal diameter, as was 

recently solved by Horvay and Poritsky [5], is another special case in which 

2 X Ag = @, Y, = Y, = dosh’, o, = |, gg = —!] 


The analytic functions are 


_— 
whe 


PG - { 


a iwR® .. 
log (¢ — 1) ae 


1 


twh> l tw” bis BP 
i ae log (¢° — 1) - 9 (wh 


The normal stress components at any point on the boundary of the disk except the two 
the supports are given by 
wR 


sin 0 


(7, ,_=9 (r = 


00) p=1 
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REPRESENTATION OF NONLINEAR FIELD FUNCTIONS BY THIELE 
SEMI-INVARIANTS' 


By WILLIAM M. MacDONALD, III,? JOHN M. RICHARDSON, ? 
AND LEON P. ROSENBERRY* 


1. Certain nonlinear field functions, possessing the property of space localization of 
the gradients of the dependent variables, occur in the problem of flame propagation in 
contmuous media. The difficulties encountered in solving the equations of propagation 
may be considerably diminished by introducing new dependent variables (intrinsically 
connected with this localization) within a certain region, outside of which the equations 
may be linearized and treated in a point-wise sense by well-known methods. Specifically, 
we choose the z-axis as nearly perpendicular to the flame front and define the region as 
x,(y, 2) S x S x.(y, 2). Taking a state variable, temperature 7’, for example, we choose 
the new dependent variables as the x, given by 


> (it)’ = = log [(¢2)] (1.1) 


Vv: 


4 


where 
i 
mc (1 2) 


Ox 


g(t) = [ de 


These definitions are closely related to the formalism of mathematical statistics. In 
particular, if one considers 07'/dx to correspond to an unnormalized distribution function 
with range (xz, S x S 22), the x, correspond to the Thiele semi-invariants [1] of 07'/dz 
and completely describe a given function in the range (7, S x S 22). 
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The semi-invariants were chosen (in preference to the moments 


[dee (1.3) 


Ox 


for example) because evaluation of the integrals 
is : oT 
, Ox 


(1.4) 


led to application of the Fourier inversion formula to the characteristic function (1.2), 
where the exponential form considerably lessened the difficulties involved. The investiga- 


tion was extended to include integrals of the form 


‘~ 7 mM m Ohms 
| de2"T,T, +++ T,—=* - 
Pas Ox 
2. The characteristic function ¢(f) of 07/dx, in the region 2,(y, z) S x S 2x2(y, 2) is 
given by Eq. (1.2). Using the Fourier inversion formula one finds 


I 


P 7 WT 
— | dte“¢)=—, mS282 (2.1) 
am . Or 

= 0), otherwise, 


from which the following expression for 7’ can be obtained if d7'/dx satisfies certain mild 


restrictions [2]: 


T(x) — T(x,) = (+) dt(e"'?’ —e''**)t '(t) dt. (2.2) 
2n/ J_. 
Writing 
oie VA 
d(te'* = a2) exp [7t(& — x,)], (2.3) 
i 0& 


it can be seen that, for 87'/dx bounded, ¢()e ‘‘” has no poles or other singularities other 
than an essential singularity at infinity and that here 
d(tje '* —~O0ast—o~, with —7/2 < argt < 7/2. 
We use this fact to reduce Eq. (2.2) to a more suitable form by writing 


ax 


| die’? —e'™)t'¢() =P | dip(te'*t' — P | dip(tje "tt", (2.4) 
where P {=., --- denotes the Cauchy proper value 
Lim | -- | | 
«0 e x ve : 


Contour integration and application of the theory of residues then yields for the last 
integral on the right of (2.4) the expression 


a x 


P dt e~"t"'6(t) = ing(0). (2.5) 
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If ¥(x) is defined as 


then from Eqs. (1. 


(2) = 


2), (2.2) and (2.: 
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T(x) — 4[7(x,) + T(2.)], (2.6) 


5), one obtains 


= (+) | dt e''*t"'¢(t) (2.7) 


}(x) 
ea 2 
and 
O34 (x) l e* it ° 
a = | dt ¢ “o(t). (2.8) 
Ox ME? dae 
In Eqs. (2.7) and (2.8) and henceforth we conventionally denote P f*. by f* 
3. We now apply (2.7) and (2.8) to 
, j ve ay 
¢,(t) = | dee" Sy — 3.1) 
. Ox 


obtaining 


“ 


o.() = ian" | d 


Using independent variables ¢, , ¢, , 


Li) 
4 


Since 


8 /¢ -s-~-l 
p(t) = 2°(2n) 


To evaluate 


one notes that 


o”’(0) 


and hence that 


| 
: ae 


rel {[ die 71 b(t) [ are “oo |, (3.2) 


_ t, , one can write 


| dt [] [t, b(t.) @(to) exp ize, , 


« 0 


ww 
iw) 
as 


where o,(t) = t — y ic « 
1 


' f2.. dx e'™ = 6(u), one can write 


[ di, «-- | dt, [] ¢'o(t,)o(c,). (3.4) 


= (d’o Ge \ean S Gy = Vite s 


¢,(t) = >, a,t"/v!. (3.7) 


v= 


The a, cannot be obtained in analytical form but one can easily expand a, in a Taylor’s 


series 


Mites , hse oes ee 


= (a,)p + , 
= 


— (da, le< 0a, ae 
( ) ~ a 9 eR ho ( ca Ki K; + ea (3.8) 
=3 P = 1=3 a 


OK; [23 OK; OK; 
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about the point P(x, x, K2, 0, 0,0 ---) in «space (x; defined by Eq. (1.1)). The coeffi- 


cients in this expansion are found to be 


(a,)p = Ab.o(—1)"(x2/2)""5,H,(7,), (3.9) 
(322) ~ Ase] 8, ot ( p 2 6)s.0} (3.10) 
OK,/ p p=1 


and 


[2 At 3] 8 — ss * > gE + &,, > ie + Sol , s) |, (3.11) 
u/ P " p=1 p= ots 


where the following notation has been introduced: 


As, = (t/2n)'t""/olu!, 


5, mi [ dt, +: [ dt, [| [o(t,) u6( > ) an operator, 
: 1 q=1 


Sy. = [{€obp(€)/dr($)}"’ Jreo , 


and //,(7,) is the vth Hermite polynomial. 

!. Treatment of the equations of hydrodynamics by the moment transformation will 
clearly involve integrals in which the integrand includes various combinations of the 
gradient of a dependent variable and some other variable. We must therefore generalize 
the foregoing treatment to the case where {% is not merely the temperature but a general- 
ized vector in a space of n dimensions; this fact can be represented symbolically by writing 
ee (4.1) 


= (& * ** 
Se a 


the product of two vectors 3 and 7 gives a vector in the product space of § and 9, 


or in oul notation. 


(4.2) 


Designating by W(t) the characteristic function of the component &; of 3, one finds, 


as In section 2, that from 
: (7t)” es sen Ota ' 
y(t) = exp { ps Ki» | = | dice"* - ; (4.3) 


follows 


03; l j atzr ‘ 
= _ (1) | | dt y,(0), (4.4) 


and 


3. = ( i) [aurea (4.5) 
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Here, %; is again the deviation from the average over the infinite range of x of the corre- 
sponding dependent variable. The canonical integral 


sy 


Fai Fa he 
| dice | ~ ( 1.6) 
, Ox 


represents then the characteristic function ¢(/) whose components are given by the 
equation 


: an 
o die (i) = ge BS i a (4.7 
2 m+. | Wa - Wam Ox ? ) 
which is easily transformed by (4.5) into the form 
f ; \ » . 0 m 
ee ft) = ( ) | dt. «+. | dt..( I] v..(t) t,) Vans) (4.8) 
277 ‘ ; , 
where £,, = ¢ — pe t, . We write, as before, 
De divcmct - aoe tee il, (4.9) 
where 


114,°** 


a, aide = Qa, ,a3,°°* we (|) 8 


the coefficients of this series then being found by a Taylor’s series expansion of the form 


about the point in x-space for which x 


0 for n = 3. The first term is given by 


Ogata De = Ws, 66,200 mai iO) fP (4.11) 
which is, from (4.8), 


er ce ‘ “ (W.,(t;)' 
\Q, Ir >= (=) | dt, | dt,, I] ( 7 7 [v2.6 | . ‘ (4.12) 


On carrying out the operations indicated on the last factor in the integrand one 
finds (see appendix | 


l 


ae)» = C, | dr( [1 4) lexp n't — YAOI), (4.13) 


. 


Oa, fie y : P 
( = ) = *] dr ti{exp (ix’t — t’/At)}H,(7,), ix<xm-+1, (4.14) 
\O Kaiu/ P Le ‘ 


Ow \ C “" ‘ [ - : . - ~ 
ae -( a J dr} I] 4] ftexp (ix’t — 'Ad}S., (4.15) 
ee a Ke: Js j=1 


and 
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where 


. m v/2 m+1 
y U vf Kam+1.2 : } 
C,= () (—1) ( 5 ) exp > Kio - 


The exponential is the matrix representation of the quadric, 


it~ ite <2 D> (kei2 + Kemer. — : ae se? 


9 
- = p>a=l 


+ i > Kaj1b; — Kan411 > “|, 


j=1 1=1 


1/2 m ° 
Kam+12 t Kam+1.1 
™ = ( 9 ) ( > > t, + ce t), 
= p=1 Kam+1.2 
and the integral is over m-dimensional spce. 


APPENDIX I—Evaluation of (a,)p and (da,/0xk,) p 
These integrals (3.9, 3.10, 3.11) all reduce to integrals of the form 


mi Ma mn 


| dr exp [—Af(x, , 22, °°* , Xn) |ar' a2" ++ Tn", 
with >> mx even if n is odd and vice versa (m; integers), 


(as 5 fay *** she) = > a+ ; bP % 2; « 
i 1 


1) 


Obviously, several integrations by parts would reduce this to an evaluation of a well- 
known integral. A simpler technique is available. For example, with s = 2 we insert the 
parameters a, 6 to form: 


a & a oO 


f(a, b) = | dx, | dx.(2,22)' exp [— (ax; + x; + ba,x2)]. 


One then evaluates f°} db af/db to find 
f(a, b) = —2x sin™ (b/2a'”’). 
Setting a = b = 1 in the derivatives of this function then gives the values of the required 


integrals. 
Moments up to the fifth have been tabulated for s = 1 and s = 2. For s = 3 only first 


terms for moments up to the fourth have been evaluated. 
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ON A ‘‘PARADOX” IN BEAM VIBRATION THEORY* 


By E. H. LEE (Brown University) 


1. Introduction. Timoshenko’ discusses the problem of a constant vertical point force 
P moving with constant velocity v across a uniform elastic beam of length / simply 
supported at its ends at the same level as shown in Figure 1. The beam is considered 
to be at rest and in equilibrium when the force commences to move across the span at 
t = 0. The force sets the beam in vertical oscillation, more or less violent depending on 
the magnitude of the force and its speed of traverse. In general the beam will be in 
oscillation, involving both kinetic and elastic strain energy, when the force reaches the 
right hand support. The “paradox” arises concerning the source of this energy, since 
the vertical force undergoes no net vertical displacement and so might be considered to 
do no net work. 

Timoshenko’s explanation consists in considering a different but related problem. 
He states that we should consider the force to act through a frictionless constraint, so 
that the force must always be directed normally to the instantaneous direction of the 
beam at the point of action. This introduces a horizontal component of force on the 
beam which is prescribed when the motion due to the vertical force has been determined. 
This horizontal component does work, and Timoshenko shows that in the case of reso- 
nance, when the traverse time is half the fundamental period of the beam, this work is 
almost equal to the energy in the fundamental mode of oscillation, the discrepancy 
being attributed to the higher modes of oscillation. 

Although this explanation seems reasonable from a purely mechanical standpoint, it 
remains puzzling as to why one can’t specify the dynamical problem in terms of a vertical 
force. The work of Lagrange showed how useful the concept of frictionless constraints 
can be in dynamics, but this does not imply that one must use them in specifying a 
problem instead of prescribing forces. The careful discussion of energy relations given 
below shows that the paradox does not in fact exist even with a prescribed vertical force. 
Since the same type of difficulty arises occasionally in other problems, it was thought 
worthwhile to write this note to clarify the issue. 

2. Energy relations. This discussion is concerned with the usual linear theory of 
small lateral vibration of beams. We shall use Timoshenko’s notation, which in addition 
to the quantities specified in Fig. 1 specifies: Young’s modulus E, second moment of 


p 


k— vf —a 


A> LO - 
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I { ’ , , 











1952 E. H. LEE 291 
the section 1, section area A, density y, and acceleration due to gravity g. 

When a force acts on a moving body, the rate of work done on the body is equal to 
the scalar product of the force and the particle velocity of the body at the point of 
application. In the present example, therefore, the rate of work done by the vertical 
is P dy/dt. However, the vertical component of the velocity of the point of 


force P 
application of the force which moves along the beam is given by a convected derivative: 
24 Oy Oy ‘ 
& = Yr yS (1 
(# a ot + Ox ) 
Since the force leaves the beam at the level at which it entered the span 
al/v / 
dy 
| (4u) dt = 0 (2) 
Jo ae 


The total work done on the beam by the constant vertical force P is: 
al/¢ >..\ 
P | (24) dt (3) 
Jo Ot] s-v1 
which in general will not be equal to zero. In fact it follows directly from (1) and (2) 
that the work (3) done by the vertical force is equal to the work done by the horizontal 
component of the constraint force in Timoshenko’s modified problem. Thus: 


i fay" nize (ay 
P | (24) dt= —P | (21) v dt. (4) 
Je OCs sues Jo ae 
The equality between the work (3) and the energy of the motion induced in the 


beam follows directly from the differential equation of the motion, and the boundary 


and initial conditions: 


,0y , yA dy , 
EI — — = P 3(z/vt) 
OX g dal 


9° 
“tu@ e288 «326 (5) 


ax” 


Oy 
y= > = 0, t= 0, 0O<r<l 
, ot es 
where 4(x/vé) is the Dirac delta function which is non zero at x = vt. Multiplying both 
sides by dy/dt and integrating with respect to x from 0 to 1 we obtain, after integration 
by parts: 


a|EI f' oy) yA f' ‘uy’ (2 ’ 
J) dy a >| = Pi ( 


Thus the rate of increase of kinetic plus strain energy is equal to the rate of work 
as defined in (3). 

3. Conclusion. Thus when energy relations are written down carefully, the fact that 
the vertical force leaves the beam at the level at which it entered the span does not 
imply no net work done. The reason for misinterpretation is that in most mechanical 
problems the point of action of an applied force does move with the body on which it 
is acting. However, there are practical examples where this is not the case. Consider 
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for example an induction motor. The magnetic force producing the torque moves with 
the synchronous angular velocity of the alternating current, whereas the rotor speed 
falls below this by the slip, so that the magnetic force is moving relative to the rotor on 
which it is producing a torque. The rate of work done on the rotor is the torque multiplied 
by the rotor speed, and not the torque multiplied by the synchronous speed. In many 
cases, as in this case, the difference in these two quantities may be lost as mechanical 
work; in this case it appears as eddy current loss. The energy balance at the point of 
application is a function of the detailed method of application of the force. A problem, 
in which a similar difficulty in the use of partial and convected derivatives arose, ap- 
peared in the discussion of the motion of a bar containing a hinge moving along its 
length.’ In this case the wrong choice involved an apparent paradox: the failure of the 
momentum principle. 


NOTE ON THE LEAST EIGENVALUE OF THE HILL EQUATION 


By TOSIO KATO (University of Tokyo) 


Let us consider the Hill equation 
y’’ + [A + f(xy = 0, —xax< x4 <om, (1) 


where f(x) is a real-valued periodic function of period 1 with the Fourier series 
f(z) ~ > c, exp (2rinz), C, = C_, » 


Recently Wintner({5], Eq. (23)) deduced the following inequalities satisfied by the lower 
limit \, of the spectrum of (1): 


—t) = = —G — 2 Dd, |e, |’. (2) 
n i 


Also the question is raised (Putnam [3], p. 314) whether the coefficient 2 on the right- 
hand side is the least possible value. In the present note we shall show that better esti- 
mates do exist. In particular, we shall show that 


] ' ‘ 
A 2 —m — = Dd Ie, |”. (3) 
For this purpose we note that \, is characterized as the least eigenvalue of (1) con- 
sidered on the finite interval 0 S x S 1 with the periodic boundary conditions 
y(0) = y(1), y’(0) = y'(1) (4) 
(see e.g. Strutt [4], p. 15). Therefore, according to the Ritz variational principle, , is 
the minimum value of the expression 


al al ; 
J[y] = (y”* — fy’) az / | y dz, 


#0 
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where y changes over real-valued functions subjected to the conditions (4). Let z(x) be 
any periodic function of period 1 with integrable z’(x). Then we have 


ol 


| (Qzeyy’ + 2’y’) dx = [zy’], = 0 


“0 


since both y and z are periodic. Hence 
| (y® —fy)dx =] [y tea? + —2 — fyy’|dz 


1 
> Min (2’ — 2 — f)- [ y dx. 


It follows that J[y] 2 Min (z’ — z° — f) and this implies 


Ayo 2 Min (2’ — z — f). (5) 
Incidentally, this is an adaptation of Wintner’s condition ({5], p. 368) to the eigenvalue 
problem under consideration. By different choices of z we can obtain different lower 


bounds of \ 
First take as z an indefinite integral g of f; = f — ¢) . Then zis certainly periodic and 


(5) becomes 
Ay = Min (—cy — g’) = —eo — (Max | g |)’. (6) 
To estimate Max | g |, we introduce the oscillation A of g: 
A = Max g — Min g. 
Since g is continuous, there are values a, b of x such that g(a) = Min g, g(b) = Max g. 
Since g is periodic, we may assume a < b < a + 1 without loss of generality. Then we 


have 
ab pat+l 
A = g(b) — g(a) = fi dz = -| fi dx. 
Ja Jb 


Hence 
As] |fild, As] |filaz. 


Addition of both inequalities and application of the Schwarz inequality give 


sa+l ol al 1/2 
2A<s | fi, |dz= | fiijdxs | fi ax | ; 


Heretofore g has been determined only up to an arbitrary additive constant. If we choose 
this constant appropriately, we can make g(b) = —g(a) = A/2. Then we have Max | g| = 
A/2 and hence 


' ' a lo, : 
(Max Pa (a= a 
‘Zig = 16 | Ji da 8 p> C. 


which, combined with (6), proves (3). 
Another estimation of \, is obtained by adjusting the arbitrary constant in g in such 
a way that 


g(x) = p (2rin) ~'c, exp (27inz). 


n#¥0 
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Then we have 


and (6) gives 


Sy ~— a ¢ ). (4) 


r = 
i 
It will be noted that (3) and (7) are mutually independent. In any case, however, (7) 
is better than the second inequality of (2), for we have 
~ ' , , , l 12 
r( La" lal sr? Da* Diel=— Dial 
) 
by the Cauchy inequality. In the same way it is seen that (7) is better than (3) if c, = 0. 


It is not clear whether the coefficient 1/8 in (3) can still be replaced by a smaller one. 
However, it cannot be made smaller than 1/2x° (Cf. Putnam [3], p. 314, where the 
figure 1/47” is given). In fact, consider the case f(x) = 2c, cos 2rx with a real c, (Mathieu 
equation); then the formula of the usual perturbation theory (Courant-Hilbert [1], p. 


300) yields easily the expansion 
Ny, = —(2z°) G + 
which is certainly convergent for sufficiently small value of | ¢, | (Kato [2], p. 169). 
It will be noted that the lower bounds of \, as given by the formulas (3) and (7) are, 
though rigorous, not very accurate from the practical standpoint. Especially this is the 


case when the Fourier coefficients c, are large, for we have 


Ao 2 —Maxf 2 —q — 2 > hs (S) 


as is easily seen by setting z = 0 in (5). 
For more accurate estimation of \, in individual cases, it is more convenient to use 
(5) directly. For instance consider the case f = 2c, cos 2rz stated above. If we set 
z = ksin 2rz, we have by 5) 
\)o = Min [(2rk — 2c,) cos 2xx — k° sin’ 27x] 
: ‘ r l 2 2 1 2 
Min [(k cos 2x7xa + er — k c,) —k (rx —k c,)"] 
k* — (wf — key)" 
This is true for every k. If we assume c, > 0 and take k = c,’", we obtain 
s 9-p'/2 2 
Ayo 2 —2e, + 2re,° — 7. 


It is easily seen that the right-hand side coincides with the asymptotic expansion of )o 
for ¢ up to the order c;”* inclusive (Strutt [4], p. 37, Eq. (4), where we have to set 
A= ry,hk = 7 C,,m, = 1). 
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